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ABSTRACT 

The Planck nominal mission cosmic microwave background (CMB) maps yield unprecedented constraints on primordial non-Gaussianity (NG). 
Using three optimal bispectrum estimators, separable template-fitting (KSW), binned, and modal, we obtain consistent values for the primordial 
local, equilateral, and orthogonal bispectrum amplitudes, quoting as our final result /^ L ral = 2.7 ± 5.8, f^ = -42 ± 75, and = -25 ± 39 
(68% CL statistical); and we find the Integrated-Sachs-Wolfe-lensing bispectrum expected in the ACDM scenario. The results are based on 
comprehensive cross-validation of these estimators on Gaussian and non-Gaussian simulations, are stable across component separation techniques, 
pass an extensive suite of tests, and are confirmed by skew-Cf, wavelet bispectrum and Minkowski functional estimators. Beyond estimates of 
individual shape amplitudes, we present model-independent, three-dimensional reconstructions of the Planck CMB bispectrum and thus derive 
constraints on early-Universe scenarios that generate primordial NG, including general single-field models of inflation, excited initial states (non- 
Bunch-Davies vacua), and directionally-dependent vector models. We provide an initial survey of scale-dependent feature and resonance models. 
These results bound both general single-field and multi-field model parameter ranges, such as the speed of sound, c s > 0.02 (95% CL), in an 
effective field theory parametrization, and the curvaton decay fraction r D > 0.15 (95% CL). The Planck data put severe pressure on ekpyrotic/cyclic 
scenarios. The amplitude of the four-point function in the local model t nl < 2800 (95% CL). Taken together, these constraints represent the highest 
precision tests to date of physical mechanisms for the origin of cosmic structure. 

Key words, cosmology: cosmic background radiation - cosmology: observations - cosmology: theory - cosmology: early Universe - cosmology: 
inflation 
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1. Introduction 

This paper, one of a set associated with the 2013 release of 
data from the Planck 1 mission (Planck Collaboration I 2013), 
describes the constraints on primordial non-Gaussianity (NG) 
obtained using the cosmic microwave background (CMB) maps 
derived from the data acquired by Planck during its nominal op- 
erations period, i.e., between 12 August 2009 and 27 November 
2010. 

Primordial NG is one of the most informative finger- 
prints of the origin of structure in the Universe, prob- 
ing physics at extremely high energy scales inaccessi- 
ble to laboratory experiments. Possible departures from 
a purely Gaussian distribution of the CMB anisotropies 
provide powerful observational access to this extreme 
physics (Allen et al. 1987; Salopek & Bond 1990; Falketal. 
1993; Gangui et al. 1994; Verde et al. 2000; Gangui & Martin 
2000b; Wang & Kamionkowski 2000; Komatsu & Spergel 
2001; Acquaviva et al. 2003; Maldacena 2003; Babich et al. 
2004; for recent reviews Bartolo et al. 2004a, Liguori et al. 
2010, Chen 2010b, Komatsu 2010, Yadav & Wandelt 2010). 
A robust detection of primordial NG - or a strong constraint 
on it - discriminates among competing mechanisms for the 
generation of the cosmological perturbations in the early 
Universe. Different inflationary models, firmly rooted in modern 
theoretical particle physics, predict different amplitudes, shapes, 
and scale dependence of NG. As a result, primordial NG is 
complementary to the scalar-spectral index of curvature pertur- 
bations and the tensor-to-scalar amplitude ratio, distinguishing 
between inflationary models that are degenerate on the basis 
of their power spectra alone. Even in the simplest models of 
inflation, consisting of a single slowly-rolling scalar field, a 
small (but calculable) level of NG is predicted (Acquaviva et al. 
2003; Maldacena 2003); this is undetectable in present-quality 
CMB and large-scale structure measurements. However, as 
demonstrated by a large body of work in recent years, extending 
this simplest paradigm will generically lead to detectable levels 
of NG in CMB anisotropies. Critically, a robust detection 
of primordial NG would rule out all canonical single-field 
slow-roll models of inflation, pointing to physics beyond the 
simplest "textbook" picture of inflation. Conversely, significant 
improvements in the constraints on primordial NG strongly 
limit extensions to the simplest paradigm, thus providing 
powerful clues to the physical mechanism that generated cosmic 
structure. 

If the primordial fluctuations are Gaussian-distributed, then 
they are completely characterised by their two-point correlation 
function, or equivalently, their power spectrum. If they are non- 
Gaussian, there is additional statistical information in the higher- 
order correlation functions, which is not captured by the two- 
point correlation function. In particular, the 3-point correlation 
function, or its Fourier counterpart, the bispectrum, is impor- 
tant because it is the lowest-order statistic that can distinguish 
between Gaussian and non-Gaussian perturbations. One of the 
main goals of this paper is to constrain the amplitude and shape 
of primordial NG using the angular bispectrum of the CMB 
anisotropies. The CMB angular bispectrum is related to the pri- 



1 Planck (http://www.esa.int/Planck) is a project of the 
European Space Agency (ESA) with instruments provided by two sci- 
entific consortia funded by ESA member states (in particular the lead 
countries France and Italy), with contributions from NASA (USA) and 
telescope reflectors provided by a collaboration between ESA and a sci- 
entific consortium led and funded by Denmark. 



mordial bispectrum defined by 

minmhWh)) = {Inf^HM+^ + k^ihMM). (1) 

Here we define the potential O in terms of the comoving cur- 
vature perturbation ( on super-horizon scales by €> = (3/5)f. 
In matter domination, on super-horizon scales, <E> is equivalent 
to Bardeen's gauge-invariant gravitational potential (Bardeen 
1980), and we adopt this notation for historical consistency. 
The bispectrum Bq,(ki,k 2 ,ki) measures the correlation among 
three perturbation modes. Assuming translational and rotational 
invariance, it depends only on the magnitudes of the three 
wavevectors. In general the bispectrum can be written as 

B< s> (k 1 ,k 2 ,k 3 ) = f NL F(kuk 2 ,k 3 ). (2) 

Here, /nl is the so-called "nonlinearity parameter" 
(Gangui et al. 1994; Wang & Kamionkowski 2000; 
Komatsu & Spergel 2001; Babich et al. 2004), a dimen- 
sionless parameter measuring the amplitude of NG. The 
bispectrum is measured by sampling triangles in Fourier space. 
The dependence of the function F(k\,k2,k-s) on the type of 
triangle (i.e., the configuration) formed by the three wavevec- 
tors describes the shape of the bispectrum (Babich et al. 
2004), which encodes much physical information. It can 
also encode the scale dependence, i.e., the running, of the 
bispectrum (Chen 2005c). 2 Different NG shapes are linked 
to distinctive physical mechanisms that can generate such 
non-Gaussian fingerprints in the early Universe. For example, 
the so-called "local" NG (Gangui et al. 1994; Verde et al. 
2000; Wang & Kamionkowski 2000; Komatsu & Spergel 2001) 
is characterized by a signal that is maximal for "squeezed" 
triangles with k\ <K k 2 — k% (or permutations; Maldacena 
2003) which occurs, in general, when the primordial NG is 
generated on super-horizon scales. Conversely, "equilateral" 
NG (Babich et al. 2004) peaks for equilateral configurations 
k\ =s k 2 ~ &3, due to correlations between fluctuation modes 
that are of comparable wavelengths, which can occur if the 
three perturbation modes mostly interact when they cross the 
horizon approximately at the same time. Other relevant shapes 
include the so-called "folded" (or flattened) NG (Chen et al. 
2007b), which is due to correlations between perturbation 
modes that are enhanced for k\ * 2k 2 m 2k-$, or the "orthogonal" 
NG (Senatore et al. 2010) that generates a signal with a positive 
peak at the equilateral configuration and a negative peak at the 
folded configuration. 

We now sketch how non-Gaussian information in the ini- 
tial conditions are transferred to observable quantities (in this 
instance, the CMB anisotropies) in the context of inflation. 
Primordial perturbations in the inflaton field(s) <p(x, t) - <f>o(t) + 
6(p(x, t) (where 5<p denotes quantum fluctuations about the 
background value 0o(O) can be characterized by the comov- 
ing curvature perturbation since this is conserved on super- 
horizon scales for adiabatic perturbations. The inflaton fluctu- 
ations 5(f) (in the flat gauge) induce a curvature perturbation 3 

2 Specifically, one can define the shape of the bispectrum as the de- 
pendence of F(ki,k 2 ,ki)(kik 2 k3) 2 on the ratios of momenta, e.g., (k 2 /ki) 
and (k^/ki), once the overall scale of the triangle K = k t + k 2 + k 3 is 
fixed. The scale dependence of the bispectrum can be characterized by 
the dependence of F(kj , k 2 , k 3 )(k] k 2 k^) 2 on the overall scale K, once the 
ratios {k 2 jk{) and (k^/ki) are fixed (see, e.g., Chen 2010b). 

3 For the curvature perturbation, we follow the notation and sign 
conventions of Komatsu et al. (2011). ( is also sometimes denoted 
ft (see e.g., Lidsey et al. 1997, Lyth & Riotto 1999 and references 
therein), while the comoving curvature perturbation ft as defined, e.g., 
in Malik & Wands (2009) is such that ft = 
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£ = -{H/(/)q)6(P at linear order; however, nonlinearities induce 
corrections to this relation. The primordial NG in the curvature 
perturbation f is intrinsically nonlinear, so that its contribution 
to the CMB anisotropics is transferred linearly at leading order. 
In particular, at the linear level, the curvature perturbation £ is 
related to Bardeen's gravitational potential O during the matter- 
dominated epoch by O = (3/5)^ and AT/T ~ g(, where g is the 
linear radiation transfer function; thus, any primordial NG will 
be transferred to the CMB even at linear order. For example, in 
the large-angular scale limit, the linear Sachs-Wolfe effect reads 
AT/T = -0/3 = -(15. Further, any other field excited during 
the inflationary phase which develops quantum fluctuations con- 
tributing to the primordial curvature perturbation - whether or 
not it is driving inflation - can leave its non-Gaussian imprint in 
the CMB anisotropies. 

Thus the bispectrum of Eq. (1) measures the fundamental 
(self-) interactions of the scalar field(s) involved in the infla- 
tionary phase and/or generating the primordial curvature per- 
turbation, as well as measuring nonlinear processes occurring 
during or immediately after inflation. It therefore brings in- 
sights into the fundamental physics behind inflation, possibly 
allowing for the first time a reconstruction of the inflationary 
Lagrangian itself. For example, in a large class of inflationary 
models which involve additional light field(s) different from the 
inflaton, the super-horizon evolution of the fluctuations in the ad- 
ditional field(s) and their transfer to the adiabatic curvature per- 
turbations can generate a large primordial NG of the local type. 
This is the case for curvaton-type models (Linde & Mukhanov 
1997; Lyth & Wands 2002; Lyth et al. 2003) where the late-time 
decay of a scalar field, belonging to the non-inflationary sec- 
tor of the theory, induces curvature perturbations; models where 
the curvature perturbation is generated by the local fluctua- 
tions of the inflaton's coupling to matter during the reheating 
phase (Kofman 2003; Dvali et al. 2004a); and multi-field models 
of inflation (see, e.g., Bartolo et al. 2002, Bernardeau & Uzan 
2002, Vernizzi & Wands 2006, Rigopoulos et al. 2006, 2007; 
Lyth & Rodriguez 2005, Byrnes & Choi 2010). Since the non- 
linear processes take place on super-horizon scales, the form of 
NG is local in real space and thus, in Fourier space, the bis- 
pectrum correlates large and small Fourier modes. "Equilateral" 
NG (Babich et al. 2004) is a generic feature of single-field mod- 
els with a non-canonical kinetic term, which can also gener- 
ate the "orthogonal" type of NG (Senatore et al. 2010). In gen- 
eral, these models are characterized by higher-derivative inter- 
actions of the inflaton field. The correlation between the fluctu- 
ation modes is suppressed when one of the modes is on super- 
horizon scales, because the derivative terms are redshifted away, 
so that the correlation is maximal for three modes of compa- 
rable wavelengths that cross the horizon at the same time. An 
example of "folded" NG is the one generated in a class of 
single-field models with non-Bunch-Davies vacuum (Chen et al. 
2007b; Holman & Tolley 2008). Indeed, these and other types of 
primordial NG can also be produced in other models, and we re- 
fer to Sect. 2 for more details. All these models can easily yield 
primordial NG with an amplitude much bigger than the one pre- 
dicted in the standard models of single-field slow-roll inflation, 
for which the NG amplitude turns out to be proportional to the 
usual slow-roll parameters / NL ~ 0(e, rj) (Acquaviva et al. 2003; 
Maldacena 2003). 

Given that a robust detection of primordial NG would 
represent a breakthrough in the understanding of the physics 
governing the Universe during its very first stages, it is crucial 
that all sources of contamination are sufficiently understood to 
firmly control their effects. In particular, any nonlinearity in 



the post-inflationary Universe can introduce NG into perturba- 
tions that were initially Gaussian. Therefore, one must ensure 
that a primordial origin is not ascribed to a non-primordial 
contaminant; however, estimators of (primordial) NG from 
CMB data will also typically be sensitive to such contaminat- 
ing signals. Potential non-primordial sources of NG can be 
classified into four broad categories: instrumental systematic 
effects (see e.g., Donzelli et al. 2009); residual foregrounds 
and point sources; secondary CMB anisotropies, such as the 
Sunyaev-Zeldovich (SZ) effect (Zeldovich & Sunyaev 1969), 
gravitational lensing (see Lewis & Challinor 2006 for a re- 
view), the Integrated Sachs-Wolfe (ISW) effect (Sachs & Wolfe 
1967) or the Rees-Sciama effect (Rees & Sciama 1968); and 
effects arising from nonlinear (second-order) perturbations 
in the Boltzmann equations (due to the nonlinear nature 
of General Relativity and the nonlinear dynamics of the 
photon-baryon fluid at recombination). Among the secondary 
anisotropies, the cross-correlation of the ISW/Rees-Sciama 
and lensing (Goldberg & Spergel 1999) produces the domi- 
nant contamination to the (local) primordial NG. The impact 
is mainly on the local type of primordial NG, because the 
ISW-lensing correlation couples the large-scale gravita- 
tional potential fluctuations sourcing the ISW effect with 
the small-scale lensing effects of the CMB, thus producing 
a bispectrum which peaks on the squeezed configurations, 
as for the local shape. Detailed analyses have shown that 
the ISW-lensing bispectrum can introduce a bias to local 
primordial NG, while the bias to equilateral primordial NG 
is negligible (see Serra & Cooray 2008, Smith & Zaldarriaga 
201 1, Hanson et al. 2009b, Lewis et al. 201 1, Mangilli & Verde 
2009, Junk & Komatsu 2012, Lewis 2012, Mangilli et al. 2013). 
In our analysis we have carefully accounted for this effect 
(we report the values of the ISW-lensing bias in Sect. 5.2, and 
demonstrate the detection of the effect with skew-Q s), as well as 
validating our results through an extensive suite of simulations 
and null tests in order to quantify the effects of systematic effects 
and diffuse and point-source foregrounds. Finally, a consistent 
treatment of weak NG in the CMB must account for additional 
contributions that arise at the nonlinear (second-order) level both 
in the gravitational perturbations after inflation ends, and for the 
evolution of the CMB anisotropies at second-order in perturba- 
tion theory at large and small angular scales. It has been shown 
that these second-order CMB effects yield negligible contami- 
nation to primordial NG for P/ancfc-quality data (Bartolo et al. 
2004b; Creminelli & Zaldarriaga 2004b; Bartolo et al. 2005; 
Boubekeur et al. 2009; Nitta et al. 2009; Senatore et al. 
2009; Khatri & Wandelt 2009; Bartolo & Riotto 2009; 
Khatri & Wandelt 2010; Bartolo et al. 2010c; Creminelli et al. 
2011b; Bartolo et al. 2012; Huang & Vernizzi 2013; Su et al. 
2012; Pettinari et al. 2013). 

Previous constraints on various shapes of primordial NG 
come from the WMAP-9 data (Bennett et al. 2012). For the lo- 
cal shape they find /^° L cal = 37 ± 20 (68% CL). For equilateral- 
type NG, they obtain /^ uil = 51 + 136 (68% CL), while for the 
orthogonal shape /°[*° = -245 ± 100 (68% CL). Other analy- 
ses employing different estimators give compatible constraints. 
Limits on other shapes, such as e.g. flattened and feature models, 
have also been obtained (Fergusson et al. 2012). 

Before concluding this section let us point out the connec- 
tion between the analyses presented here and in the compan- 
ion paper Planck Collaboration XXIII (2013) on the statistical 
and isotropy properties of the CMB. Statistical anisotropy and 
NG are essentially two alternative descriptions of the same phe- 
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nomenon on the sky (Ferreira & Magueijo 1997). Specifically 
any Gaussian but statistically anisotropic model becomes, af- 
ter averaging over the possible (a priori unknown) orienta- 
tions of the anisotropy, a statistically isotropic non-Gaussian 
model. For example local NG can be generated by large-scale 
field fluctuations that couple to the small-scale power. For 
the given fixed realization of large-scale modes that we see, 
the small-scale anisotropies look anisotropic on the sky, and 
it is equally valid to describe this as a Gaussian anisotropic 
model (assuming the large-scale modes are Gaussian). In this 
paper we mostly focus on the non-Gaussian interpretation of 
various physically motivated models, although it is useful to 
bear both perspectives in mind, in particular when considering 
what forms of non-primordial signal might cause contamina- 
tion. Planck Collaboration XXIII (2013) consider a broad class 
of more general phenomenological forms of anisotropy, which 
are complementary to the analysis presented here. 

This paper is organized as follows. In Sect. 2, we present 
models generating primordial NG that have been tested in this 
paper. Section 3 summarizes the statistical estimators used to 
constrain the CMB bispectrum from Planck data and the meth- 
ods for the reconstruction of the CMB bispectrum. Section 4 
summarizes the statistical estimator used to constrain the CMB 
trispectrum. In Sect. 5, we discuss the non-primordial contribu- 
tions to the CMB bispectrum and trispectrum, including fore- 
ground residuals after component separation and focusing on the 
/ml bias induced by the ISW-lensing bispectrum. Section 6 de- 
scribes an extensive suite of tests performed on realistic simula- 
tions to validate the different estimator pipelines, and compare 
their performance. Using simulations, we also quantify the im- 
pact on /nl of using a variety of component-separation tech- 
niques. Section 7 contains our main results: we present con- 
straints on /nl for the local, equilateral, and orthogonal bispec- 
tra, and a selected set of other bispectrum shapes; we show a 
reconstruction of the CMB bispectrum, and give limits on the 
CMB trispectrum. In Sect. 8 we validate these results by per- 
forming a series of null tests on the data to assess the robustness 
of our results. In Sect. 9, we discuss the main implications of 
Planck's constraints on primordial NG for early Universe mod- 
els. We conclude in Sect. 10. Appendix A contains a deriva- 
tion of the expected scatter between /nl results on the same 
map from different estimators used in the validation tests of 
Sect. 6, while Appendix B presents a comparison of constraints 
on some selected non-standard bispectrum shapes using different 
foreground-cleaned maps. 

2. Inflationary models for primordial 
non-Gaussianity 

There is a simple reason why standard single-field models of 
slow-roll inflation predict a tiny level of NG, of the order of the 
usual slow-roll parameters / NL ~ 0(e, r/): 4 in order to achieve 

4 This has been shown in the pioneering research which demon- 
strated that perturbations produced in single-field models of slow-roll 
inflation are characterized by a low-amplitude NG (Salopek & Bond 
1990; Falketal. 1993; Gangui et al. 1994). Later Acquaviva et al. 
(2003) and Maldacena (2003) obtained a complete quantitative pre- 
diction for the nonlinearity parameter in single-field slow-roll infla- 
tion models, also showing that the predicted NG is characterized 
by a shape dependence which is more complex than suggested by 
previous results expressed in terms of the simple parameterization 
<D(x) = <S> L (x) + /nl^M (Gangui et al. 1994; Verde et al. 2000; 
Wang & Kamionkowski 2000; Komatsu & Spergel 2001) ,where <I> L is 
the linear gravitational potential. 



an accelerated period of expansion, the inflaton potential must 
be very flat, thus suppressing the inflaton (self-)interactions and 
any sources of nonlinearity, and leaving only its weak gravita- 
tional interactions as the main source of NG. This fact leads to 
a clear distinction between the simplest models of inflation, and 
scenarios where a significant amplitude of NG can be generated 
(e.g., Komatsu 2010), as follows. The simplest inflationary mod- 
els are based on a set of minimal conditions: (i) a single weakly- 
coupled neutral single scalar field (the inflaton, which drives 
inflation and generates the curvature perturbations); (ii) with a 
canonical kinetic term; (iii) slowly rolling down its (featureless) 
potential; (iv) initially lying in a Bunch-Davies (ground) vacuum 
state. In the last few years, an important theoretical realization 
has taken place: a detectable amplitude of NG with specific tri- 
angular configurations (corresponding broadly to well-motivated 
classes of physical models) can be generated if any one of the 
above conditions is violated (Bartolo et al. 2004a; Liguori et al. 
2010; Chen 2010b; Komatsu 2010; Yadav & Wandelt 2010): 

- "local" NG, where the signal peaks in "squeezed" triangles 
(k\ <s &2 — £3) (e.g., multi-field models of inflation); 

- "equilateral" NG, peaking for k\ a* £2 ~ £3- 
Examples of this class include single-field models with 
non-canonical kinetic term (Chen et al. 2007b), such as k- 
inflation (Armendariz-Picon et al. 1999; Chen et al. 2007b) 
or Dirac -Born-Infield (DBI) inflation (Silverstein & Tong 
2004; Alishahiha et al. 2004), models characterized by more 
general higher-derivative interactions of the inflaton field, 
such as ghost inflation (Arkani-Hamed et al. 2004), and 
models arising from effective field theories (Cheung et al. 
2008); 

- "folded" (or flattened) NG. Examples of this class in- 
clude: single-field models with non-Bunch-Davies vac- 
uum (Chen et al. 2007b; Holman & Tolley 2008) and models 
with general higher-derivative interactions (Senatore et al. 
2010; Bartolo et al. 2010a); 

- "orthogonal" NG which is generated, e.g., in single- 
field models of inflation with a non-canonical kinetic 
term (Senatore et al. 2010), or with general higher-derivative 
interactions. 

All these models naturally predict values of |/nlI » 1. A de- 
tection of such a signal would rule out the simplest models of 
single-field inflation, which, obeying all the conditions above, 
are characterized by weak gravitational interactions with | /nlI ^ 
1. 

The above scheme provides a general classification of infla- 
tionary models in terms of the corresponding NG shapes, which 
we adopt for the data analysis presented in this paper: 

1. "general" single-field inflationary models (tested using the 
equilateral, orthogonal and folded shapes); 

2. multi-field models of inflation (tested using the local shape). 

In each class, there exist specific realizations of inflationary 
models which are characterized by the same underlying phys- 
ical mechanism, generating a specific NG shape. We will inves- 
tigate these classes of inflationary models by constraining the 
corresponding NG content, focusing on amplitudes and shapes. 
We also perform a survey of non-standard models giving rise 
to alternative specific shapes of NG. Different NG shapes are 
observationally distinguishable if their cross-correlation is suf- 
ficiently low; almost all of the shapes analysed in this paper 
are highly orthogonal to each other (e.g., Babich et al. 2004; 
Fergusson & Shellard 2007). 
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There are exceptional cases which evade this classification: 
for example, some exotic non-local single-field theories of in- 
flation produce local NG (Barnaby & Cline 2008), while some 
multi-field models can produce equilateral NG, e.g., if some 
particle production mechanism is present (examples include 
trapped inflation Green et al. 2009, and some models of axion 
inflation Barnaby & Peloso 2011; Barnaby et al. 2011, 2012b). 
Another example arises in a class of multi-field models where 
the second scalar field is not light, but has a mass m « H, of 
the order of the Hubble rate during inflation. Then NG with 
an intermediate shape, interpolating between local and equi- 
lateral, can be produced - "quasi-single field" models of in- 
flation (Chen & Wang 2010a,b) - for which the NG shape is 
similar to the so-called constant NG of Fergusson & Shellard 
(2007). Furthermore, there is the possibility of a superposition 
of shapes (and/or running of NG), generated if different mech- 
anisms sourcing NG act simultaneously during the inflation- 
ary evolution. For example, in multi-field DBI inflation, equi- 
lateral NG is generated at horizon crossing from the higher- 
derivative interactions of the scalar fields, and it adds to the 
local NG arising from the super-horizon nonlinear evolution 
(e.g., Langlois et al. 2008a,b; Arroja et al. 2008; Renaux-Petel 
2009). 

In the following subsections, we discuss each of these possi- 
bilities in turn. The reader already familiar with this background 
material may skip to Sect. 3. 



where Pq,(k) — A/k . 8 is the power spectrum of Bardeen's grav- 
itational potential with normalization A 2 and scalar spectral in- 
dex « s . For example, the models introduced in the string theory 
framework based on the DBI action (Silverstein & Tong 2004; 
Alishahiha et al. 2004) can be described within the P{X, </>)-class, 
and they give rise to an equilateral NG with an overall amplitude 
Jnl ~ -(35/108)c,7 2 for c, « 1, which turns out typically to 
beC u "<-5. 5 

The equilateral shape emerges also in models characterized 
by more general higher-derivative interactions, such as ghost in- 
flation (Arkani-Hamed et al. 2004) or models within effective 
field theories of inflation (Cheung et al. 2008; Senatore et al. 
2010; Bartolo et al. 2010a). 

Taken individually, each higher-derivative interaction of the 
inflaton field generically gives rise to a bispectrum with a 
shape which is similar - but not identical to - the equilateral 
form (an example is provided by the two interaction terms dis- 
cussed above for an inflaton with a non-standard kinetic term). 
Therefore it has been shown, using an effective field theory ap- 
proach to inflationary perturbations, that it is possible to build a 
combination of the corresponding similar equilateral shapes to 
generate a bispectrum that is orthogonal to the equilateral one, 
the so-called "orthogonal" shape. This can be approximated by 
the template (Senatore et al. 2010) 



2.1. General single-field models of inflation 



Typically in models with a non-standard kinetic term (or more 
general higher-derivative interactions), inflaton perturbations 
propagate with an effective sound speed c s which can be smaller 
than the speed of light, and this results in a contribution to 
the NG amplitude /nl ~ cj 2 in the limit c s <K 1. For exam- 
ple, models with a non-standard kinetic term are described by 
an inflaton Lagrangian £, - P(X, <p), where X - g^'d^ d v <p, 
with at most one derivative on <p, and the sound speed is c 2 = 
(dP/dX)/(dP/dX + 2X(d 2 P/dX 2 )). 

In this case, two interaction terms give the dominant con- 
tribution to primordial NG, one of the type (<5</>) 3 and the other 
of the type 6(p(V5<p) 2 , which arise from expanding the P(X, 0) 
Lagrangian. Each of these two interaction terms generates a 
bispectrum with a shape similar to the equilateral type, with 
the first inflaton interaction yielding a nonlinearity parameter 
/nl ~ c^ 2 , independent of the amplitude of the other bispec- 
trum. Equilateral NG is usually generated by derivative interac- 
tions of the inflaton field; derivative terms are suppressed when 
one perturbation mode is frozen on super-horizon scales during 
inflation, and the other two are still crossing the horizon, so that 
the correlation between the three perturbation modes will be sup- 
pressed, while it is maximal when all the three modes cross the 
horizon at the same time, which happens for k\ m k 2 ~ k 3 . 

The equilateral type NG is well approximated by the tem- 
plate (Creminelli et al. 2006) 



B^(k u k 2 M) = (>A 2 fZ* 
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(4) 



The orthogonal bispectrum can also arise as the predomi- 
nant shape in some inflationary realizations of Galileon infla- 
tion (Renaux-Petel et al. 201 1). 



Non-separable single-field bispectrum shapes: While most 
single-field inflation bispectra can be well-characterized by the 
equilateral and orthogonal shapes, we note that these are sep- 
arable ansatze which only approximate the contributions from 
two leading order terms in the cubic Lagrangian. In an effective 
field theory approach these correspond to two shapes which can 
be associated directly with the inflaton field interactions 7r(<9,7r) 2 
and 7r 3 (in the language of the effective field theory of inflation 
the inflaton scalar degree of freedom n is related to the comov- 
ing curvature perturbation as ( — -Hn). They are, respectively 



5 An effectively single-field model with a non-standard kinetic term 
and a reduced sound speed for the adiabatic perturbation modes might 
also arise in coupled multi-field systems, where the heavy fields 
are integrated out: see discussions in, e.g., Tolley & Wyman (2010); 
Achucarro et al. (2011); Shiu & Xu (2011). 
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(Senatore et al. 2010, see also Chen et al. 2007b) 
B™ T \k u k 2 M) 



tEftI/i. , ,,_6A 2 / N E r 1 (-9/17) 



Z*f+ZNV^- 3 *^1 
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(6) 



These shapes differ from equilateral in the flattened or collinear 
limit. DBI inflation gives a closely related shape of particular 
interest phenomenologically (Alishahiha et al. 2004), 



->DB! , / , , 6A2 C' (-3/7) 



B^(k u k 2 ,k 3 ) 



(kik 2 k 3 ) 3 {k x + k 2 + k 3 ) 2 
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(7) 



For brevity, we have given the scale-invariant form of the shape 
functions, without the mild power spectrum running. There are 
also sub-leading order terms which give rise to additional non- 
separable shapes, but these are expected to be much smaller 
without special fine-tuning. 

2.2. Multi-field models 

This class of models generally includes an additional light 
scalar field (or more fields) during inflation, which can be 
different from the inflaton, and whose fluctuations contribute to 
the final primordial curvature perturbation of the gravitational 
potential. It could be the case of inflation driven by several 
scalar fields - "multiple-field inflation" - or the one where the 
inflaton drives the accelerated expansion, while other scalar 
fields remain subdominant during inflation. This encompasses, 
for instance, a large class of multi-field models which leads to 
non-Gaussian isocurvature perturbations (for earlier works, see 
e.g., Linde & Mukhanov 1997, Peebles 1997, Bucher & Zhu 
1997). More importantly, such models can also lead to cross- 
correlated and non-Gaussian adiabatic and isocurvature modes, 
where NG is first generated by large nonlinearities in some 
scalar (possibly non-inflatonic) sector of the theory, and 
then efficiently transferred to the inflaton adiabatic sector(s) 
through the cross-correlation of adiabatic and isocurvature 
perturbations 6 (Bartolo et al. 2002; Bernardeau & Uzan 
2002; Vernizzi & Wands 2006; Rigopoulos et al. 2006, 
2007; Lyth & Rodriguez 2005; Tzavara & van Tent 2011; 
for a review on NG from multiple-field inflation models, 
see, Byrnes & Choi 2010). Another interesting possibility 
is the curvaton model (Mollerach 1990; Enqvist & Sloth 
2002; Lyth & Wands 2002; Moroi & Takahashi 2001), where 



6 This may happen, for instance, if the inflaton field is coupled 
to the other scalar degrees of freedom, as expected on particle 
physics grounds. These scalar degrees of freedom may have large self- 
interactions, so that their quantum fluctuations are intrinsically non- 
Gaussian, because, unlike the inflaton case, the self-interaction strength 
in such an extra scalar sector does not suffer from the usual slow-roll 
conditions. 



a second light scalar field, subdominant during inflation, 
decays after inflation ends, producing the primordial den- 
sity perturbations which can be characterized by a high NG 
level (e.g., Lyth & Wands 2002; Lyth et al. 2003; Bartolo et al. 
2004d). NG in the curvature perturbation can be generated 
at the end of inflation, e.g., due to the nonlinear dynamics of 
(p)reheating (e.g., Enqvist et al. 2005; Chambers & Rajantie 
2008; Barnaby & Cline 2006; see also Bond et al. 2009) or, as 
in modulated (p)reheating and modulated hybrid inflation, due 
to local fluctuations in the decay rate/interactions of the inflaton 
field (Kofman 2003; Dvali et al. 2004a,b; Bernardeau et al. 
2004; Zaldarriaga 2004; Lyth 2005; Salem 2005; Lyth & Riotto 
2006; Kolb et al. 2006; Cicoli et al. 2012). The common feature 
of all these models is that a large NG in the curvature pertur- 
bation can be produced via both a transfer of super-horizon 
non-Gaussian isocurvature perturbations in the second field (not 
necessarily the inflaton) to the adiabatic density perturbations, 
and via additional nonlinearities in the transfer mechanism. 
Since, typically, this process occurs on super-horizon scales, 
the form of NG is local in real space. Being local in real 
space, the bispectrum correlates large and small scale Fourier 
modes. The local bispectrum is given by (Falketal. 1993; 
Gangui et al. 1994; Verde et al. 2000; Wang & Kamionkowski 
2000; Komatsu & Spergel 2001) 

B^\k x ,k 2 ,k 3 ) = 2/i L cal [p (fc 1 ) J P fe) + P (^ 1 )P (A:3) 
+ P«(fc)Po(Jfc 3 )] 

1 



_ ry a 2 /-local 

- ZA Jnl 



,4-n s .4-« s 



+ cycl. 



(8) 



Most of the signal-to-noise ratio in fact peaks in the squeezed 
configurations (k\ <K k 2 £3) 



0,k 2 ,k 3 ) 



(9) 



The typical example of a curvature perturbation that generates 
the bispectrum of Eq. (8) is the standard local form for the 
gravitational potential (Hodges et al. 1990; Kofman et al. 1991; 
Salopek&Bond 1990; Gangui et al. 1994; Verde et al. 2000; 
Wang & Kamionkowski 2000; Komatsu & Spergel 2001) 



<D(x) 



<D L (x) + f^^lix) - <<!>£(*)» , 



(10) 



where Q>l(x) is the linear Gaussian gravitational potential and 
f^ 1 is the amplitude of a quadratic nonlinear correction (though 
this is not the only possibility: e.g., the gravitational potential 
produced in multiple-field inflation models generally cannot be 
reduced to the Eq. (10)). For example, in the (simplest) adiabatic 
curvaton models, the NG amplitude turns out to be (Bartolo et al. 
2004d,c) f§f = (5/4r D ) - 5r D /6 - 5/3, for a quadratic po- 
tential of the curvaton field (Lyth & Wands 2002; Lyth et al. 
2003; Lyth & Rodriguez 2005; Sasaki et al. 2006), where r D = 
[3pcurvaton/(3p C urvaton + 4p rad iation)]D is the "curvaton decay frac- 
tion" evaluated at the epoch of the curvaton decay in the sudden 
decay approximation. Therefore, for ru <s 1, a high level of NG 
is imprinted. 

There exists a clear distinction between multi-field and 
single-field models of inflation that can be probed via a con- 
sistency condition (Maldacena 2003; Creminelli & Zaldarriaga 
2004a; Chen et al. 2007b; Chen 2010b): in the squeezed limit, 
single-field models predict a bispectrum 



^single-field, 7 
B a, (M 



0,k 2 ,k 3 =*2)' 



-(l-nMkMM) , (11) 



6 



Planck Collaboration: Planck 2013 Results. XXIV. Constraints on primordial NG 



and thus /nl ~ 0(n s - 1) in the squeezed limit, in a model- 
independent sense (i.e., not only for standard single-field mod- 
els). This means that a significant detection of local NG (in the 
squeezed limit) would rule out a very large class of single-field 
models of inflation (not just the simplest ones). Although based 
on very general conditions, the consistency condition of Eq. (11) 
can be violated in some well-motivated inflationary settings (we 
refer the reader to Chen (2010b); Chen et al. (2013) and refer- 
ences therein for more details). 

Quasi-single field inflation: Quasi-single field inflation has an 
extra field (or fields) with mass m close to the Hubble parame- 
ter H during inflation; these models evolve quiescently, produc- 
ing a calculable non-Gaussian signature (Chen & Wang 2010b). 
The resulting one-parameter bispectrum smoothly interpolates 
between local and equilateral models, though in a non-trivial 



Bf\k u k 2 ,h) 



3^ 2 N v [8k l k 2 k i /(k l + k 2 + k 3 ) 3 \ 
(hk 2 k 3 yi 2 N v [8/27](A:i +k 2 + fc 3 ) 3/2 



(12) 



where v = (9/4 - m 2 /H 2 ) 1 ^ 2 and N v is the Neumann function 
of order v. Quasi-single field models can also produce an es- 
sentially "constant" bispectrum defined by B const (k\, k 2 , k 3 ) = 
6A 2 /^™' st I '(k\k 2 ki) 2 . The constant model is the simplest possible 
non-zero primordial shape, with all its late-time CMB structure 
simply reflecting the behaviour of the transfer functions. 

Alternatives to inflation: Local NG can also be generated 
in some alternative scenarios to inflation, for instance in 
cyclic/ekpyrotic models (for a review, see Lehners 2010), due 
to the same basic curvaton mechanism described above. In this 
case, typical values of the nonlinearity parameter can easily 
reach > 10. 

2.3. Non-standard models giving rise to alternative specific 
forms of NG 

Non-Bunch-Davies vacuum and higher-derivative interactions: 
Another interesting bispectrum shape is the folded one, which 
peaks in flattened configurations. To facilitate data analyses, 
the flat shape has been usually parametrized by the tem- 
plate (Meerburg et al. 2009) 

BZ(k u k 2 ,k 3 ) = 6A 2 f^ 



,4-« s ,4-n s 



,4-« s ,4-« s 
K 2 K 3 



k\~ n *k\~ n * 
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(13) 



The initial quantum state of the inflaton is usually specified 
by requiring that, at asymptotically early times and short dis- 
tances, its fluctuations behave as in flat space. Deviations from 
this standard "Bunch-Davies" vacuum can result in interesting 
features in the bispectrum. Models with an initial non-Bunch- 
Davies vacuum state (Chen et al. 2007b; Holman & Tolley 
2008; Meerburg et al. 2009) can generate sizeable NG similar 
to this type. NG highly correlated with such a template can 
be produced in single-field models of inflation from higher- 
derivative interactions (Bartolo et al. 2010a), and in models 
where a "Galilean" symmetry is imposed (Creminelli et al. 
2011a). In both cases, cubic inflaton interactions with two 



derivatives of the inflaton field arise. Single-field inflation 
models with a small sound speed, studied in Senatore et al. 
(2010), can generate the flat shape, as a result of a linear 
combination of the orthogonal and equilateral shapes. In fact, 
from a simple parametrization point of view, the flat shape 
can be always written as i*fl a t(£i, £2, £3) = [Feq U ii(k\,k 2 ,k 3 ) - 
F nho(k\,k 2 ,kj,)]/2 (Senatore et al. 2010). Despite this, we pro- 
vide constraints also on the amplitude of the flat bispectrum 
shape of Eq. (13). 

For models with excited (i.e., non-Bunch-Davies) initial 
states, the resulting NG shapes are model-dependent, but they 
are usually characterized by the importance of flattened or 
collinear triangles, with £3 » k\ + k 2 along the edges of the 
tetrapyd. We will denote the original flattened bispectrum shape, 
given in Eq. (3.62) of Chen et al. (2007b), by B™ BD ; it is gener- 
ically much more flattened than the "flat" model of Eq. (13). 
Although this shape was derived specifically for power-law k- 
inflation, it encapsulates several different shapes, with ampli- 
tudes which can vary between different phenomenological mod- 
els. These shapes are also typically oscillatory, being regular- 
ized by a cutoff scale k c giving the oscillation period; this cutoff 
k c m (c s T c y l is determined by the (finite) time t c in the past 
when the non-Bunch-Davies component was initially excited. 
For excited canonical single-field inflation, the two leading order 
shapes can be described (Agullo & Parker 201 1) by the ansatz 
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(14) 



where f\{k\,k 2 ,k 3 ) = k 2 (k 2 + k 2 )/2 is dominated by squeezed 
configurations, f 2 (k\,k 2 , k 3 ) = k 2 k 2 has a flattened shape, and i = 
1 , 2. Note that for all oscillatory shapes, the relevant bispectrum 
equation defines the normalisation of /nl- The flattened signal 
is most easily enhanced in the limit of small sound speed c s , for 
which a regularized ansatz is given by (Chen et al. 2007b) 
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NBD3 
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(15) 



Scale-dependent feature and resonant models: Oscillating bis- 
pectra can be generated from violation of a smooth slow-roll 
evolution ("feature" or "resonant" NG). These models have the 
distinctive property of a strong running NG, which breaks ap- 
proximate scale-invariance. A sharp feature in the inflaton po- 
tential forces the inflaton field away from the attractor solu- 
tion, and causes oscillations as it relaxes back; these oscillations 
can appear in the bispectrum (Wang & Kamionkowski 2000; 
Chen et al. 2007a, 2008), as well as the power spectrum and 
other correlators. An analytic form for the oscillatory bispectrum 
for these feature models is (Chen et al. 2007a) 



B f ^{h,k 2 ,h) = 



6A 2 / r 



NL 



(hk 2 hf 



2n{k\ + k 2 + £3) 



3k c 



(16) 



where <p is a phase factor and k c is a scale associated with the 
feature, which is linked in turn to an effective multipole period- 
icity i c of the CMB bispectrum. Typically, these oscillations will 
decay with an envelope of the form exp[-(fei + k 2 + k 3 )/mk c ] for 
a model-dependent parameter m. 

Closely related "resonant" bispectra can be created by pe- 
riodic features superimposed on a smooth inflation potential 
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(Chen et al. 2008; Flauger & Pajer 2011); these induce small 
periodic features in the background evolution, with which the 
quantum inflaton fluctuations can resonate while still inside 
the horizon. Resonant models are particularly relevant in the 
context of axion inflation models (e.g., Flauger et al. 2010; 
Flauger & Pajer 2011; Barnaby et al. 2012b). These mecha- 
nisms also create oscillatory behaviour in the bispectrum, but 
with a more constant amplitude and a wavelength that becomes 
logarithmically stretched. Here, the resonant oscillations for 
most models can be represented in the form 

6A 2 F 68 

B™(k u k 2 ,k 3 ) = - ■ NL . sin [Cln(h +k 2 + k 3 ) + <f>] , (17) 
(kik 2 k 3 y 

where the constant C = 1/ ln(3A; c ) and <p is a phase. 

Finally, we note that periodic features in the inflationary po- 
tential can excite the vacuum state, as well as perturbing the 
background inflation trajectory (Chen 2010a). Such models offer 
the intriguing possibility of combining the flattened non-Bunch- 
Davies shape with periodic oscillations: 



B™" BD (k u k 2 ,k 3 ) 



t a 2 rresNBD 



{exp(-£ c 3/5 (£ 2 +£3-£i)/2fci) 



x sin[£ c ((fc 2 + h - k\)llk\ + lnfci) + <p] + 2 perm.} . (18) 

This ansatz represents the dominant folded resonant contribution 
in inflationary models with non-canonical kinetic terms, which 
competes with resonant (Eq. ( 1 7)) and equilateral (Eq. (3)) con- 
tributions; however, for slow-roll single-field inflation, there are 
additional terms. 

Directional dependence motivated by gauge fields: Additional 
variations of the bispectrum shape have been proposed for mod- 
els with vector fields, which can have an additional direc- 
tional dependence through the parameter p l2 - k\ ■ k 2 where 
k = k/k. For example, primordial magnetic fields sourcing 
curvature perturbations can cause a dependence on both p and 
p 2 (Shiraishi et al. 2012), and a coupling between the inflaton 
<p and the gauge field strength F 2 can yield a p 2 dependence 
(Barnaby et al. 2012a; Bartolo et al. 2013). We can parameter- 
ize these shapes as variations on the local shape, following 
Shiraishi et al. (2013), as 

B^ikuhM = Yj c dPLUin)Pa>{h)P<s>{k 2 ) + 2 perm], (19) 

L 

where Pl(p) is tne Legendre polynomial with Pq = 1, P\ = p. 
and P 2 - ^(p 2 - 1). For example, for L — 1 we have the shape 
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Also the recently introduced "solid inflation" 
model (Endlich et al. 2012) generates bispectra similar to 
Eq. (19). Here and in the following the nonlinearity parameters 
f£ h are related to the c L coefficients by cq = 2f^°,c\ - — 



16/^ L 2 . The L = 1,2 shapes exhibit sharp variations 



NL 

and c\ 

in the flattened limit for e.g., k\ + k 2 * k 3 , while in the squeezed 
limit, L = 1 is suppressed whereas L — 2 grows like the local 
bispectrum shape (i.e., the L — case). Whether or not the 
underlying gauge field models prove robust, this directional 
dependence on the wave vectors is a generic feature which 
yields distinct bispectrum families, deserving closer study. 

Warm inflation: In warm inflation (Berera 1995), where dissipa- 
tive effects are important, a non-Gaussian signal can be gener- 
ated (e.g., Moss & Xiong 2007) that peaks in the squeezed limit 



- but with a more complex shape than the local one - and ex- 
hibiting a low cross-correlation with the other shapes (see refer- 
ences in Liguori et al. 2010). 

2.4. Higher-order non-Gaussianity: the trispectrum 

The connected four-point functions of CMB anisotropics (or the 
harmonic counterpart, the so-called trispectrum) can also pro- 
vide crucial information about the mechanism that gave rise to 
the primordial curvature perturbations (Okamoto & Hu 2002). 
The primordial trispectrum is usually characterised by two am- 
plitudes tnl and g^L- tnl * s most often related to / 2 L -type con- 
tributions, while <7nl is the amplitude of intrinsic cubic nonlin- 
earities in the primordial gravitational potential (corresponding, 
in terms of field interactions, to a scalar-exchange and to a con- 
tact interaction term, respectively). They correspond to 'soft' 
limits of the full four-point function, with respectively the di- 
agonal and one side of the general wavevector trapezoid being 
much smaller than the others. In the CMB maps they appear re- 
spectively approximately as a spatial variation in amplitude of 
the small-scale fluctuations, and a spatial variations in the value 
of /nl correlated with the large-scale temperature. In addition to 
possible primordial signals that are the focus of this paper there 
is also expected to be a large lensing trispectrum (of very dif- 
ferent shape), discussed in detail in Planck Collaboration XVII 
(2013). 

The simplest local trispectrum is given by 

<®(Jfci)a>(Jfc 2 )0(Jk 3 )*(*4)> = (2n) 3 6 < - 3 \k l +k 2 + k 3 + k 4 ) 
(25 

x I — T NL [P<b(ki)P<s,(k 2 )Ps,(kn) + (11 perm.)] 



+6<?nl [PMPMP^ki) + (3 perm.)] 



(21) 



where fey = \k t + kj\. Previous constraints on r NL and gr N LS 
have been derived, e.g., by Smidt et al. (2010) who obtained 
-7.4 x 10 5 < g NL < 8.2 x 10 s and -0.6 x 10 4 < t nl < 3.3 x 10 4 
(at 95% CL) analysing WMAP-5 data; for the same datasets 
Fergusson et al. (2010b) obtained -5.4 x 10 5 < g Nh < 8.6 X 10 5 
(68% CL). This kind of trispectrum typically arises in multi-field 
inflationary models where large NG arise from the conversion of 
isocurvature perturbations on superhorizon scales. If the curva- 
ture perturbation is the standard local form, in real space one has 
0(x) = d> L (x) + f^W L (x) - <02» + g NL (S>l(x). In this case, 
tnl = /5) 2 ; however, in general the trispectrum ampli- 

tude can be larger. 

The trispectrum is a complementary observable to the CMB 
bispectrum as it can further distinguish different inflationary sce- 
narios. This is because the same interactions that lead to the bis- 
pectrum might be responsible also for a large trispectrum, so 
that the different NG parameters can be related to each other in 
a well-defined way within specific models. If there is a non-zero 
squeezed-shape bispectrum there must necessarily be a trispec- 
trum, with T NL > (6/^ L cal /5) 2 (Suyama & Yamaguchi 2008; 
Sugiyama et al. 201 1; Sugiyama 2012; Lewis 201 1; Smith et al. 
2011; Assassi et al. 2012; Kehagias & Riotto 2012). In the sim- 
plest inflationary scenarios the prediction would be tnl = 
{6f^ al /5) 2 , but larger values would indicate more complicated 
dynamics. Several inflationary scenarios have been found in 
which the bispectrum is suppressed, thus leaving the trispec- 
trum as the largest higher-order correlator in the data. A detec- 
tion of a large trispectrum and a negligible bispectrum would 
be a smoking gun for these models. This is the case, for ex- 
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ample, of certain curvaton and multi-field field models of in- 
flation (Byrnes et al. 2006; Sasaki et al. 2006; Byrnes & Choi 
2010), which for particular parameter choice can produce a 
significant tnl and 0nl and small /nl- Large trispectra are 
also possible in single-field models of inflation with higher- 
derivative interactions (see, e.g., Chen et al. 2009; Arroja et al. 
2009; Senatore & Zaldarriaga 2011; Bartolo et al. 2010b), but 
these would be suppressed in the squeezed limit since they 
are generated by derivative interactions at horizon-crossing, and 
hence only project weakly onto the local shapes. These equi- 
lateral trispectra arise can be well-described by some template 
forms (Fergusson et al. 2010b). Naturally, higher-order correla- 
tions could also be considered, but are not directly studied in this 
paper. 



3. Statistical estimation of the CMB bispectrum 

In this Section, we review the statistical techniques that we use 
to estimate the nonlinearity parameter / NL . We begin by fixing 
some notation and describing the CMB angular bispectrum in 
Sect. 3.1. We then introduce in Sect. 3.2 the optimal /nl bispec- 
trum estimator. From Sect. 3.2.1 onwards we describe in detail 
the different implementations of the optimal estimator that were 
developed and applied to Planck data. 



3. 1 . The CMB angular bispectrum 

Temperature anisotropics are represented using the a.f m coeffi- 
cients of a spherical harmonic decomposition of the CMB map, 



AT sr~t 

— («) = 2_ u a tmYt m {n); 

Cm 



(22) 



we write C( = (\ac m \ 2 ) for the angular power spectrum and Q = 
(2t + \y l Yum \ a em\ 2 for the corresponding (ideal) estimator; hats 
denote estimated quantities. The CMB angular bispectrum 
is the three-point correlator of the a[ m : 



r^nim 2 m 3 _ , > 

t\£ 2 £ 3 — \ ii E\m\ u l 2 m 2 Lt l 3 m 3 / ■ 



(23) 



If the CMB sky is rotationally invariant, the angular bispectrum 
can be factorized as follows: 



{a( imi a{ 2 m 2 ae 3 m } } - Q^tmlm, 



(24) 



where £>f,f 2 f 3 is the so called reduced bispectrum, and @mimim 3 ls 
the Gaunt integral, defined as: 



_ » / (\ (2 h 

- l\hh I mi mi mj 

where he l e 2 e } is a geometrical factor. 



(2A + 1)(2€ 2 + 1)(2€ 3 + 1) / U h (3 



An 



o o • 



(25) 



(26) 



The Wigner-3 j symbol in parentheses enforces rotational sym- 
metry, and allows us to restrict attention to a tetrahedral domain 
of multipole triplets {^1,^2,^3}, satisfying both a triangle condi- 
tion and a limit given by some maximum resolution ( max (the 
latter being defined by the finite angular resolution of the ex- 
periment under study). This three-dimensional domain TV of 



allowed multipoles, sometimes referred to in the following as a 
"tetrapyd", is illustrated in Fig. 1 and it is explicitly defined by 

Triangle condition: 6\ < £2 + ^3 f° r t\ ^ ^2< ^3,+perms., 
Parity condition: €\ + h + I3 = 2n , n e N , (27) 
Resolution: t x , € 2 , h < 4™ , £\,h,&3 e N. 

Here, TV is the isotropic subset of the full space of bispectra, 
denoted by T 7 . 

One can also define an alternative rotationally-invariant re- 
duced bispectrum B[ l f 2 c i in the following way: 



mi miwjn * ' 



t\ £2 £3 J j^m\m 2 m 3 (^8) 
nil m 2 m 3 } ' _c 



Note that this Bf,M 3 is equal to hf^t, times the angle-averaged 
bispectrum as defined in the literature. From Eqs. (24) and (25), 
and the fact that the sum over all ra, of the Wigner-3 j symbol 
squared is equal to 1, it is easy to see that Z?f,f,f 3 is related to the 
reduced bispectrum by 



B 



(29) 



The interest in this bispectrum B^t^ is that it can be estimated 
directly from maximally-filtered maps of the data: 



B 



(30) 



h = J d 2 hT h (h)T (2 (h)T h {h) , 

where the filtered maps T((h) are defined as: 

Te{n) = a (m Y tm {n) . (31) 

m 

This can be seen by replacing the B^T 3 in Eq. (28) by its 
estimate at imi a{ im a^ m3 and then using Eq. (25) to rewrite the 
Wigner symbol in terms of a Gaunt integral, which in its turn 
is expressed as an integral over the product of three spherical 
harmonics. 

3.2. CMB bispectrum estimators 

The full bispectrum for a high-resolution map cannot be eval- 
uated explicitly because of the sheer number of operations in- 
volved, 0(^, 5 nax ), as well as the fact that the signal will be too 
weak to measure in individual multipoles with any significance. 
Instead, we essentially use a least-squares fit to compare the 
bispectrum of the observed CMB multipoles with a particular 
theoretical bispectrum b^^. We then extract an overall "am- 
plitude parameter" /nl for that specific template, after defin- 
ing a suitable normalization convention so that we can write 
bt&h ~ fiwJ^if , where b^ { ( is defined as the value of the 
theoretical bispectrum ansatz for / NL = 1. 

Optimal 3-point estimators, introduced by Heavens (1998) 
(see also Gangui & Martin 2000a), are those which saturate the 
Cramer-Rao bound. Taking into account the fact that instrument 
noise and masking can break rotational invariance, it has been 
shown that the general optimal /nl estimator can be written 
as (Babich 2005; Creminelli et al. 2006; Senatore et al. 2010; 
Verde et al. 2013): 



ti,mi 



fix 



(32) 



X \£tim u t\m\ a t\™\ ^ ' e 2 m 2 ,t' 2 m' 2 a l' 2 m' 2 ( ^ 't 3 m 3 ,{' 3 mf 3 a t' 3 m' 3 



3C 



-l 



C 



-l 



l\tn\,€ 2 m 2 ' l 3 m 3 ,t\m'S^£'-i m \ 
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L t\ 

Fig. 1. Permitted observational domain of Eq. (27) for the CMB bispec- 
trum b f] f 2 { 3 . Allowed multipole values £ 2 , (3) lie inside the shaded 
"tetrapyd" region (tetrahedron+pyramid), satisfying both the triangle 
condition and the experimental resolution £ < L = £ max . 



where C -1 is the inverse of the covariance matrix Cf, mi /, m , = 
(ae imi ae 2mi ) and N isa suitable normalization chosen to produce 
unit response to ( . 

In the expression of the optimal estimator above we note the 
presence of two contributions, one (hereafter defined the "cubic 
term" of the estimator) is cubic in the observed ai m and corre- 
lates the bispectrum of the data to the theoretical fitting template 

e , while the other is linear in the observed ai, n (hereafter, 
the "linear term"), which is zero on average. In the rotationally- 
invariant case the linear term is proportional to the monopole in 
the map, which has been set to zero, so in this case the estima- 
tor simply reduces to the cubic term. However, when rotational 
invariance is broken by realistic experimental features such as 
a Galactic mask or an anisotropic noise distribution, the linear 
term has an important effect on the estimator variance. In this 
case, the coupling between different t would in fact produce a 
spurious increase in the error bars (coupling of Fourier modes 
due to statistical anisotropy can be "misinterpreted" by the esti- 
mator as NG). The linear term correlates the observed ac m to the 
power spectrum anisotropies and removes this effect, thus restor- 
ing optimality (Creminelli et al. 2006; Yadav et al. 2008, 2007). 

The actual problem with Eq. (32) is that its direct imple- 
mentation to get an optimal /nl estimator would require mea- 
surement of all the bispectrum configurations from the data. As 
already mentioned at the beginning of this section, the compu- 
tational cost of this would scale like ^ ax and be totally pro- 
hibitive for high-resolution CMB experiments. Even taking into 
account the constraints imposed by isotropy, the number of mul- 
tipole triples {(1,(2,(3) is of the order of 10 9 at Planck resolu- 
tion, and the number of different observed bispectrum configu- 
rations ^mfmjU i s °f me order of 10 15 . For each of them, costly 
numerical evaluation of the Wigner symbol is also required. This 
is completely out of reach of existing supercomputers. It is then 
necessary to find numerical solutions that circumvent this prob- 
lem and in the following subsections we will show how the dif- 
ferent estimators used for the /nl Planck data analysis address 



this challenge. Before entering into a more accurate description 
of these different methods, we would like however to stress again 
that they are all going to be different implementations of the opti- 
mal /nl estimator defined by Eq. (32); therefore they are concep- 
tually equivalent and expected to produce /nl results that are in 
very tight agreement. This will later on allow for stringent vali- 
dation tests based on comparing different pipelines. On the other 
hand, it will soon become clear that the different approaches that 
we are going to discuss also open up a range of additional ap- 
plications beyond simple /nl estimation for standard bispectra. 
Such applications include, for example, full bispectrum recon- 
struction (in a suitably smoothed domain), tests of directional 
dependence of /nl, and other ways to reduce the amount of data, 
going beyond simple single-number /nl estimation. So different 
methods will also provide a vast range of complementary infor- 
mation. 

Another important preliminary point, to notice before dis- 
cussing different techniques, is that none of the estimators in 
the following sections implement exactly Eq. (32), but a slightly 
modified version of it. In Eq. (32) the CMB multipoles always 
appear weighted by the inverse of the full covariance matrix. 
Inverse covariance filtering of CMB data at the high angular 
resolutions achieved by experiments like WMAP and Planck is 
another very challenging numerical issue, which was fully ad- 
dressed only recently (Smith et al. 2009; Komatsu et al. 2011; 
Eisner & Wandelt 2013). For our analyses we developed two in- 
dependent inverse-covariance filtering pipelines. The former is 
based on an extension to Planck resolution of the algorithm used 
for WMAP analysis (Smith et al. 2009; Komatsu et al. 201 1); the 
latter is based on the algorithm described in Eisner & Wandelt 
(2013). However, detailed comparisons interestingly showed 
that our estimators perform equally well (i.e., they saturate the 
Cramer-Rao bound) if we approximate the covariance matrix as 
diagonal in the filtering procedure and we apply a simple diffu- 
sive inpainting procedure to the masked areas of the input CMB 
maps. A more detailed description of our inpainting and Wiener 
filtering algorithms can be found in Sect. 3.3. 

In the diagonal covariance approximation, the minimum 
variance estimator is obtained by making the replacement 
(C~ l a)( m — > acmlCt in the cubic term and then including the 
linear term that minimizes the variance for this class of cubic 
estimator (Creminelli et al. 2006). This procedure leads to the 
following expression: 



f — _ \ r 1 1 '2 Ci T, th v 
L ~ N Zj^'i" 1 * AW 

— ^— 6 — - — , (55) 

. Q, C( 2 Cf 3 Q, Cf 2 C(^ m 

where the tilde denotes the modification of C( and b^^j, to in- 
corporate instrument beam and noise effects, and indicates that 
the multipoles are obtained from a map that was masked and pre- 
processed through the inpainting procedure detailed in Sect. 3.3. 
This means that 

hi( 2 f 3 = b^be.b^be^ , Q = b 2 e C e + N t , (34) 

where be denotes the experimental beam, and N( is the noise 
power spectrum. For simplicity of notation, in the following we 
will drop the tilde and always assume that beam, noise and in- 
painting effects are properly included. 
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Using Eqs. (28) and (29) we can rewrite Eq. (33) in terms of 
the bispectrum Bf,f 2 f 3 : 



NL 



6_ 

N 



E 



Dth / Dobs _ nlin A 



(35) 



t 2 h 



In the above expression, B th is the theoretical template for B 
(with /ml = 1) and Z? obs denotes the observed bispectrum (the cu- 
bic term), extracted from the (inpainted) data using Eq. (30). B hn 
is the linear correction, also computed using Eq. (30) by replac- 
ing two of the filtered temperature maps by simulated Gaussian 
ones and averaging over a large number of them (three permuta- 
tions). The variance V in the inverse-variance weights is given by 
Vfj^fj = gt^t^ f ( C^,Cf 2 C< 3 (remember that these should be 
viewed as being the quantities with tildes, having beam and noise 
effects included) with ge^e, a permutation factor (ge^e, = 6 
when all i are equal, gc^( 3 - 2 when two £ are equal, and 
<fcil 2 l3 = 1 otherwise). Both Eqs. (33) and (35) will be used in the 
following. Eq. (33) will provide the starting point for the KSW, 
skew-Q and modal estimators, while Eq. (35) will be the basis 
for the binned and wavelets estimators. 

Next, we will describe in detail the different methods, and 
show how they address the numerical challenge posed by the 
necessity to evaluate a huge number of bispectrum configura- 
tions. To summarize loosely: the KSW estimator, the skew-Q 
approach and the separable modal methodology achieve mas- 
sive reductions in computational costs by exploiting structural 
properties of b lh , e.g., separability. On the other hand, the binned 
bispectrum and the wavelet approaches achieve computational 
gains by data compression of B obs . 

3.2.1. The KSW estimator 

To understand the rationale behind the KSW es- 
timator (Komatsu et al. 2005, 2003; Senatore et al. 
2010; Creminelli et al. 2006; Yadav et al. 2008, 2007; 
Yadav & Wandelt 2008; Smith & Zaldarriaga 2011), as- 
sume that the theoretical reduced bispectrum £>* w , can be 
exactly decomposed into a separable structure, e.g., there exist 
some sequences of functions a{£, r), B((, r) such that we can 
approximate fe/,^ 3 as 



?W 3 



r)B(€ 2 , r)a(h, r) + 8(€ u r)B{€ z , r)a{€ 2 , r) 



fj3(4r»,r)^i,r)] r 2 dr, 



(36) 



where r is a radial coordinate. This assumption is fulfilled 
in particular by the local shape (Komatsu & Spergel 2001; 
Babich et al. 2004), with a((, r) and B{(, r) involving integrals of 
products of spherical Bessel functions and CMB radiation trans- 
fer functions. Let us consider the optimal estimator of Eq. (33) 
and neglect for the moment the linear part. Exploiting Eq. (36) 
and the factorizability property of the Gaunt integral (Eq. (25)), 
the cubic term of the estimator can be written as: 



where 



and 



S cub = j drr 2 j d 2 hA(h, r)B 2 (h, r) 
A(n,r) = > , 

B(n ' r) = Z, c, ■ 



(37) 
(38) 
(39) 



> lin 



From the formulae above we see that the overall triple inte- 
gral over all the configurations l\, 1%, {3 has been factorized into 
a product of three separate sums over different t. This pro- 
duces a massive reduction in computational time, as the prob- 
lem now scales like ^ ax instead of the original ^ ax . Moreover, 
the bispectrum can be evaluated in terms of a cubic statistic 
in pixel space from Eq. (37), and the functions A(h, r), B(h, r) 
are obtained from the observed ci( m by means of Fast Harmonic 
Transforms. 

It is easy to see that the linear term can be factorized in anal- 
ogous fashion. Again considering the local shape type of decom- 
position of Eq. (36), it is possible to find: 

= Y J drr 2 j d 2 h [2{A(r,h)B(r,h)) MC X 
x B(r, h) + (B(r, h)B(r, h)) MC A(r, «)] , (40) 

where (-)mc denotes a Monte Carlo (MC) average over simula- 
tions accurately reproducing the properties of the actual data set 
(basically we are taking a MC approach to estimate the prod- 
uct between the theoretical bispectrum and the ac m covariance 
matrix appearing in the linear term expression). 

The estimator can be finally expressed as a function of S cu b 
and S Mn : 

/nl= ^ + 5^_ (41) 

Whenever it can be applied, the KSW approach makes the prob- 
lem of /nl estimation computationally feasible, even at the high 
angular resolution of the Planck satellite. One important caveat 
is that factorizability of the shape, which is the starting point of 
the method, is not a general property of theoretical bispectrum 
templates. Strictly speaking, only the local shape is manifestly 
separable. However, a large class of inflationary models can be 
extremely well approximated by separable equilateral and or- 
thogonal templates (Babich et al. 2004; Creminelli et al. 2006; 
Senatore et al. 2010). The specific expressions of cubic and lin- 
ear terms are of course template-dependent, but as long as the 
template itself is separable their structure is analogous to the ex- 
ample shown in this Section, i.e., they can be written as pixel 
space integrals of cubic products of suitably-filtered CMB maps 
(involving MC approximations of the a( m covariance for the lin- 
ear term). For a complete and compact summary of KSW im- 
plementations for local, equilateral and orthogonal bispectra see 
Komatsu et al. (2009, Appendix). 

3.2.2. The Skew-Q Extension 

The skew-Q statistics were introduced by Munshi & Heavens 
(2010) to address an issue with estimators such as KSW which 
reduce the map to an estimator of /nl for a given type of NG. 
This level of data compression, to a single number, has the dis- 
advantage that it does not allow verification that a NG signal is 
of the type which has been estimated. KSW on its own cannot 
tell if a measurement of /nl of given type is actually caused by 
NG of that type, or by contamination from some other source or 
sources. The skew-Q statistics perform a less radical data com- 
pression than KSW (to a function of €), and thus retain enough 
information to distinguish different NG signals. The desire to 
find a statistic which is able to fulfil this role, but which is still 
optimal, drives one to a statistic which is closely related to KSW, 
and indeed reduces to it when the scale-dependent information 
is not used. A further advantage of the skew-Q is that it allows 
joint estimation of the level of many types of NG simultaneously. 



11 



Planck Collaboration: Planck 2013 Results. XXIV. Constraints on primordial NG 



This requires a very large number of simulations for accurate es- 
timation of their covariance matrix, and they are not used in this 
role in this paper. However, they do play an important part in 
identifying which sources of NG are clearly detected in the data, 
and which are not. 

We define the skew-CV statistics by extending from KSW, as 
follows: from Eq. (37), the numerator <5 can be rewritten as 



(42) 



where 

cy 2 : 



and 



P{l x \r)a txm Y txmx (ti) 
Ct, 



P(€ 2 ; r)a tim2 Y iim2 (h) a({; r)a tm Y tm {h) 



c. 



r 2 d 2 hdr (43) 



c r - /JT.Z Z 

<U *Jd- e c_ „,, ... „ 



(ft) 



a({ 2 ; r)at 2mi Yi 2Ml (h) j3{£; r) 



iY{ m (n) 



Ct 



Ct 



r 2 d 2 hdr. (44) 



The skew-Q approach allows for the full implementation of the 
KSW procedure, when the sum in Eq. (42) is fully evaluated; 
furthermore, it allows for extra degrees of flexibility, e.g., by re- 
stricting the sum to subsets of the multipole space, which may 
highlight specific features of the NG signal. Each form of NG 
considered has its own a,/3, hence its own set of skew-CV, de- 

a d2 AB B 

noted S c = + 2Cp, and we have chosen to illustrate 

here with the local form, but as with KSW the method can be 
extended to other separable shapes, and some skew-Q do not 
involve integrals, such as the ISW-lensing skew statistic. Note 
that in this paper we do not fit the S( directly, but instead we 
estimate the NG using KSW, and then verify (or not) the nature 
of the NG by comparing the skew-Q with the theoretical expec- 
tation. No further free parameters are introduced at this stage. 
This procedure allows investigation of KSW detections of NG 
of a given type, assessing whether or not they are actually due to 
NG of that type. 

3.2.3. Separable Modal Methodology 

Primordial bispectra need not be manifestly separable (like the 
local bispectrum), or be easily approximated by separable ad 
hoc templates (equilateral and orthogonal), so the direct KSW 
approach above cannot be applied generically (nor to late-time 
bispectra). However, we can employ a highly-efficient gener- 
alization by considering a complete basis of separable modes 
describing any late-time bispectrum (see Fergusson & Shellard 
2007; Fergusson et al. 2010a), as applied to WMAP-7 data for 
a wide variety of separable and non-separable bispectrum mod- 
els (Fergusson et al. 2012). See also Planck Collaboration XXIII 
(2013) and Planck Collaboration XXV (2013). We can achieve 
this by expanding the signal-to-noise-weighted bispectrum as 



VQ1Q2Q3 u,k 



(45) 



where the (non-orthogonal) separable modes Q„ are defined by 



cyclic perms in i, j,k] . 



(46) 



It is more efficient to label the basis as Q„, with the subscript n 
representing an ordering of the [i, j, k] products (e.g., by distance 
i 2 + j 2 + k 2 ). The q t (€) are any complete basis functions up to 
a given resolution of interest and they can be augmented with 
other special functions adapted to target particular bispectra. The 
modal coefficients in the bispectrum of Eq. (45) are given by the 
inner product of the weighted bispectrum with each mode 



n ' 



p < ^C[,Cf 2 Ce } 
where the modal transformation matrix is 



(47) 



Jnp = (Qn, Qp) 

• 2> ! 



i, 1 it. 



(48) 



In the following, the specific basis functions q~j{f) we employ in- 
clude either weighted Legendre-like polynomials or trigonomet- 
ric functions. These are combined with the Sachs-Wolfe local 
shape and the separable ISW shape in order to obtain high cor- 
relations to all known bispectrum shapes (usually in excess of 
99%). 

Substituting the separable mode expansion of Eq. (45) into 
the estimator and exploiting the separability of the Gaunt integral 
(Eq. (25)), yields 

6 = ^2 Yf* n J d2fl [M\p(h)M r (h)M s] (h) 

~6(^{h)M?{K))M s) {h)\ . (49) 



n<r^prs 



where the M p (h) represent versions of the original CMB map 
filtered with the basis function q p (and the weights ( VQ)" 1 ), 
that is, 

M p (h) = q p {€) ^ YUh) • (50) 



tm 

The maps M p (h) incorporate the same mask and a realistic 
model of the inhomogeneous instrument noise; a large ensem- 
ble of these maps, calculated from Gaussian simulations, is used 
in the averaged linear term in the estimator of Eq. (49), allow- 
ing for the subtraction of these important effects. Defining the 
integral over these convolved product maps as cubic and linear 
terms respectively, 

A, cub = Jd 2 hM lp (h)M r (h)M s} (h), (51) 

P„ ]in = Jd 2 h(M°(h)M?(h))M s] (h), (52) 

the estimator reduces to a simple sum over the mode coefficients 



(53) 



where ft® = ^f t cah -ff^ m . The estimator sum in Eq. (53) is now 
straightforward to evaluate because of separability, since it has 
been reduced to a product of three sums over the observational 
maps (Eq. (49)), followed by a single 2D integral over all di- 
rections (Eq. (51)). The number of operations in evaluating the 
estimator sum is only 0{£ 2 ). 

For the purposes of testing a wide range inflationary mod- 
els, we can also define a set of primordial basis functions 
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Qi„k(ki,k,2,k-)) = cji(ki)qj(k2)qk(kx) + perms, for wavenumbers 
satisfying the triangle condition (again we will order the [i, j, k) 
with n). This provides a separable expansion for an arbitrary 
primordial shape function S(ki,k2,ks) = B(k\,k2,k2)l{k\k2k-s) 2 , 
that is, 

S (h , k 2 , k 3 ) = ^ QnQnih , k 2 , h) . (54) 

n 

Using the same transfer functions as in the KSW inte- 
gral (37), we can efficiently project forward each separa- 
ble primordial mode Q n {ki,k2,k 3 ) to a corresponding late- 
time solution Q„(l\l2h) (essentially the projected CMB bispec- 
trum for Q„(k\, k2, £3)). By finding the inner product between 
these projected modes Q„{i\,(2^i) an d the CMB basis func- 
tions Qn(A>^2!^3)> we can obtain the transformation matrix 
(Fergusson et al. 2010a,b) 

Tn P =Jr^l.«!.<3). Q(dJ2,m, (55) 
r 

which projects the primordial expansion coefficients af t to late- 
time: 

a„ = 2^ r np or p . (56) 
p 

When T np is calculated once we can efficiently convert any given 
primordial bispectrum B{k\,k2,k-i) into its late-time CMB bis- 
pectrum counterpart using Eq. (45). We can use this to extend 
the KSW methodology and to search for the much broader range 
of primordial models beyond local, equilateral and orthogonal, 
having validated on these standard shapes. 

3.2.4. Binned Bispectrum 

In the binned bispectrum approach (Bucheretal. 2010), the 
computational gains are achieved by data compression of the ob- 
served B. This is quite feasible, because like the power spectrum 
the bispectrum is a rather smooth function, with features on the 
scale of the acoustic peaks. In this way one obtains an enormous 
reduction of the computational resources needed at the cost of 
only a tiny increase in variance (typically 1%). 

More precisely, the following statistic is considered, 

+t 

(e&i m=-e 

where the A, are intervals (bins) of multipole values [i h l M - 1], 
for 2=0,..., (TYbins - 1), with t = { min and £ Nbbs = £ max + 1, and 
the other bin boundaries chosen in such a way as to minimize the 
variance of /nl- The binned bispectrum is then obtained by using 
Tj instead of 7> in the expression for the sample bispectrum of 
Eq. (30), to obtain: 

B W 3 = f T h (h)T i2 (h)T h (h)d 2 h. (58) 

The linear term B ]ln is binned in an analogous way, and the the- 
oretical bispectrum template B th and variance V are also binned 
by summing them over the values of t inside the bin. Finally 
/nl is determined using the binned version of Eq. (35), i.e., by 
replacing all quantities by their binned equivalent and replac- 
ing the sum over i by a sum over bin indices ;. An important 
point is that the binned bispectrum estimator does not mix the 
observed bispectrum and the theoretical template weights until 
the very last step of the computation of /nl (the final sum over 



bin indices). Thus, the (binned) bispectrum of the map is also 
a direct output of the code. Moreover, one can easily study the 
^-dependence of the results by omitting bins from this final sum. 

The full binned bispectrum allows one to explore the bis- 
pectral properties of maps independent of a theoretical model. 
Binned bispectra have been used to compare component separa- 
tion maps and single-frequency maps dominated by foregrounds, 
as presented in Sects. 3.4.2 and 7. In its simplest implementa- 
tion, which is used in this paper, the binned estimator uses top- 
hat filters in harmonic space, which makes the Gaussian noise 
independent between different bins; however, slightly overlap- 
ping bins could be used to provide better localization properties 
in pixel space. In this sense the binned estimator is related to the 
wavelet estimators, which we discuss below. 

3.2.5. Wavelet / N l estimator 

Wavelet methods are well-established in the CMB litera- 
ture and have been applied to virtually all areas of the 
data analysis pipeline, including map-making and compo- 
nent separation, point source detection, search for anoma- 
lies and anisotropies, cross-correlation with large scale struc- 
ture and the ISW effect, and many others (see for instance 
Antoine & Vandergheynst 1998, Martinez-Gonzalez et al. 2002, 
Cayon et al. 2003, McEwen et al. 2007a, Pietrobon et al. 2006, 
Starcketal. 2006, McEwen et al. 2007b, Cruz et al. 2007, 
Fayetal. 2008, Feeney et al. 2011, Geller & Mayeli 2009a, 
Geller & Mayeli 2009b, Starck et al. 2010, Scodeller et al. 201 1, 
Fernandez-Cobos et al. 2012). These methods have the advan- 
tage of possessing localization properties both in real and 
harmonic space, allowing the effects of masked regions and 
anisotropic noise to be dealt with efficiently. 

In terms of the current discussion, wavelets can be viewed 
as a way to compress the sample bispectrum vector by a 
careful binning scheme in the harmonic domain. See also 
Planck Collaboration XXIII (2013). In particular, the wavelet 
bispectrum can be rewritten as 

1 1 r , 

q ijk - d hw(Rj,h)w(Ri,n)w(R k ,h), (59) 

4tt (r ;( r ; ,n J 

where 

w(R; b)= f d 2 fi f(nmh- b,R) = Y, a fm co e (R)Y em (h). (60) 

Here b is the position on the sky at which the wavelet coefficient 
is evaluated and cr is is the dispersion of the wavelet coefficient 
map w(R; b) for the scale R. The filters cl>((R) can be seen as 
the coefficients of the expansion into spherical harmonic of the 
Spherical Mexican Hat Wavelet (SMHW) filter 



*P(;c,n;fl) = 




Here N(R) = R yj\ +R 2 /2 + R 4 /4 is a normalizing constant and 
y = 2 tan(0/2) represents the distance between x and n, evaluated 
on the stereographic projection on the tangent plane at «; 6 is 
the corresponding angular distance, evaluated on the spherical 
surface. 

The implementation of the linear-term correction can pro- 
ceed in analogy with the earlier cases. However, note that, in 
view of the real-space localization properties of the wavelet 
filters, the linear term here is smaller than for KSW and re- 
lated approaches, although not negligible. Moreover, it can be 
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well-approximated by a term-by-term sample-mean subtrac- 
tion for the wavelet coefficients, which allows for a further re- 
duction of computational costs. Further details can be found 
in Curtoetal. (2011a,b, 2012); Regan et al. (2013) (see also 
Lan & Marinucci 2008; Rudjord et al. 2009; Pietrobon et al. 
2009, 2010; Donzelli et al. 2012 for related needlet-based pro- 
cedures). 

3.3. Wiener filtering 

As discussed in Sect. 3.2, the /nl bispectrum estimator requires, 
in principle, inverse covariance filtering of the data to achieve 
optimality (equivalent to Wiener filtering up to a trivial multipli- 
cation by the inverse of the signal power spectrum). 

We have used the iterative method of Eisner & Wandelt 
(2013) for Wiener filtering simulations and data. The algorithm 
makes use of a messenger field, introduced to mediate between 
the pixel space and harmonic space representation, where noise 
and signal properties can be specified most directly. Since the 
Wiener filter is the maximum a posteriori solution, we moni- 
tor the x 2 of the current solution as a convergence diagnostics. 
We terminate the algorithm as soon as the improvement in the 
posterior between two consecutive steps has dropped below a 
threshold of A^ 2 < lO^oy, where cr^ is the standard deviation 
of ^-distributed variables with a number of degrees of freedom 
given by the number of observed pixels. Results of this method 
have been cross-validated using an independent conjugate gra- 
dient inversion algorithm with multi-grid preconditioning, orig- 
inally developed for the analysis of WMAP data in Smith et al. 
(2009). Applying this estimator to simulations pre-processed us- 
ing the above mentioned algorithms yielded optimal error bars 
as expected. 

However, we found that optimal error bars could also be 
achieved for all shapes using a much simpler diffusive inpaint- 
ing pre-filtering procedure that can be described as follows: all 
masked areas of the sky (both Galactic and point sources) are 
filled in with an iterative scheme. Each pixel in the mask is filled 
with the average of all surrounding pixels, and this is repeated 
2000 times over all masked pixels (we checked on simulations 
that convergence of all / NL estimates was achieved with 2000 it- 
erations). Note that the effect of this inpainting procedure, espe- 
cially visible for the Galactic mask, is effectively to apodize the 
mask, reproducing small-scale structure near the edges and only 
large-scale modes in the interior. This helps to prevent propa- 
gating any sharp-edge effects or lack of large-scale power in the 
interior of the mask to the unmasked regions during harmonic 
transforms. 

Any bias and/or excess variance arising from the inpainting 
procedure were assessed through MC validation (see Sect. 6) and 
found to be negligible. Since the inpainting procedure is partic- 
ularly simple to implement, easily allows inclusion of realistic 
correlated-noise models in the simulations, and retains optimal- 
ity, we chose inpainting as our data filtering procedure for the 
/nl analysis. 

3.4. CMB bispectrum reconstruction 

3.4.1. Modal bispectrum reconstruction 

Modal and related estimators can be used to reconstruct the 
full bispectrum from the modal coefficients y8y* obtained from 
map filtering with the separable basis functions 2yifc(A> (2^3) 
(Eq. (51)) (Fergusson et al. 2010a). It is easy to show that the 
expectation value for (or equivalently /3„) for an ensemble 



of non-Gaussian maps generated with a CMB bispectrum shape 
given by a„ (Eq. (47)) (with amplitude /nl) is 



</?„> = /nl^ 



(62) 



where y np is defined in Eq. (48). Hence, we expect the best fit 
coefficients for a particular a„ realization to be given by the /?„s 
themselves (for a sufficiently large signal). Assuming that we 
can extract the /3 n coefficients accurately from a particular ex- 
periment, we can directly reconstruct the CMB bispectrum using 
the expansion of Eq. (46), that is, 



n.p 



(63) 



where, for convenience, we have also defined orthonormalized 
basis functions R n {t\, I2, (3) with coefficients af t and such that 



<R», R P ) 



>np 



This method has been validated for simulated 



maps, showing the accurate recovery of CMB bispectra, and it 
has been applied to the WMAP-1 data to reconstruct the full 3D 
CMB bispectrum (Fergusson et al. 2012). 

To quantify whether or not there is a model-independent 
deviation from Gaussianity, we can consider the total inte- 
grated bispectrum. By summing over all multipoles, we can de- 
fine a total integrated nonlinearity parameter F NL which, with 
the orthonormal modal decomposition of Eq. (63) becomes 
(Fergusson et al. 2012) 



NL 



1 H 1 
1 n t\l 



2 b 2 



<R2 



y a R1 

Lun LI 11 



(64) 



where N\ oc is the normalization for the local /nl = 1 model. 
The expectation value (^ NL ) contains more than just the three- 
point correlator, with contributions from products of the two- 
point correlators and even higher-order contributions, 



NL 



1 



"max 

6« ma x + J] ( F NL«f + <K 2 >6p< + -) 



(65) 



Here F 2 L is the full 6-point function over the unconnected 
Gaussian part, i.e., the product of 3 Q. So, for a Gaussian in- 
put this recovers an average of 1 per mode. In the non-Gaussian 
case the leading-order contributions after the Gaussian part are 
the bispectrum squared and the product of the power spectrum 
and the trispectrum, which enter at the same order (for additional 
explanations see Fergusson et al. 2012). 



3.4.2. Smoothed binned bispectrum reconstruction 

As explained in Sect. 3.2.4, the full binned bispectrum of the 
maps under study is one of the products of the binned estimator 
code. 

Given the relatively fine binning (about 50 bins up to I - 
2000 or 2500), most of the measurement in any single bin-triplet 
is noise. If combinations of maps are chosen so that the CMB 
primordial signal dominates, most of this noise is Gaussian, re- 
flecting the fact that even when (xyz) = 0, for a particular statis- 
tical realization, xyz is almost certainly non-zero. If our goal is 
to test whether there is any statistically significant signal, then it 
makes sense to normalize by defining a new field Sjjj 2 ,- 3 , which is 
the binned bispectrum divided by its expected standard deviation 
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(computed in the standard way assuming Gaussian statistics). 
The distribution of S,,,-,,, in any one bin-triplet is very nearly 
Gaussian as a result of the central limit theorem, and there is 
almost no correlation between different bins. 

We could study the significance of the extreme values at this 
fine resolution, but it also makes sense to smooth in order to 
detect features coherent over a wider range of I. If the three- 
dimensional domain over which the bispectrum is defined (con- 
sisting of those triplets in the range [{ m \ n , Anax] satisfying the tri- 
angle inequality and parity condition) did not have boundaries, 
the smoothing would be more straightforward. We smooth with 
a Gaussian kernel of varying width in A(ln €), normalized so that 
the smoothed function in each pixel would be normal-distributed 
if the input map were Gaussian. In this way, based on the extreme 
values, it is possible to decide on whether there is a statistically 
significant NG signal in the map in a "blind" (non-parametric) 
way. 



4. Statistical estimation of the CMB trispectrum 

4.1. The squeezed-diagonal trispectrum: r NL 

The 4-point function, or equivalently the trispectrum, of the 
CMB, can also place interesting constraints on inflation- 
ary physics. There are several physically interesting "shapes" 
of the trispectrum (e.g., Huang & Shiu 2006; Byrnes et al. 
2006; Fergusson et al. 2010b; Izumi et al. 2012), in analogy 
to the bispectrum case. In the simplest non-Gaussian mod- 
els, the CMB bispectrum has larger signal-to-noise than the 
trispectrum, but there are examples of technically natural 
models in which the trispectrum has larger signal-to-noise 
(e.g., Senatore & Zaldarriaga 2011; Baumann & Green 2012; 
see also Bartolo et al. 2010b). This can happen in models in 
which the field modulating the fluctuation amplitude is only 
weakly correlated to the observed large-scale curvature pertur- 
bation. 

Analysis of the trispectrum is more challenging than that of 
the bispectrum, due to the increased range of systematic effects 
and secondary signals which can contribute. For example, grav- 
itational lensing of the CMB generates a many-sigma contribu- 
tion to the trispectrum, though it has a distinctive anisotropic 
shape that differs from primordial NG modulated by scalar 
fields. As an instrumental example, any mismatch between the 
true covariance of the observed CMB plus noise and the co- 
variance which is assumed in the analysis (due, for example, to 
mischaracterisation of the pointing, beams, or noise properties) 
will generally lead to biases in the estimated trispectrum. Due to 
these challenges, we have deferred a full analysis of the primor- 
dial trispectrum to a future paper, and here focus on the simplest 
squeezed shape that can provide useful constraints on primordial 
models, tnl- 

Tnl is most easily understood as measuring the large-scale 
modulation of small-scale power. The constraints on /nl show 
that such a modulation must be small if correlated with the tem- 
perature. However it is possible for multi-field inflation models 
to produce squeezed-shape modulations which are uncorrelated 
with the large-scale curvature perturbations. Such models can be 
constrained by the trispectrum, conventionally parameterized by 
tnl in the squeezed-diagonal shape. 

For example, consider the case where a small-scale Gaussian 
curvature perturbation £q is modulated by another field <p so that 
the primordial perturbation is given by 



where <p(x) is a large-scale modulating field (with amplitude 
<k 1). The large-scale modes of can be measured from the 
modulation they induce in the small-scale £ power spectrum. If 
<p has a nearly scale-invariant spectrum, the nearly-white cos- 
mic variance noise on the reconstruction dominates on small- 
scales, so only the very largest modes can be reconstructed 
(Kogo & Komatsu 2006). A reconstruction of <p is going to be 
limited to only very large-scale variations, in which case the 
scale of the variation is very large compared to the width of the 
last-scattering surface; i.e., in any particular direction a large- 
scale modulating field will modulate all small-scale perturba- 
tions through the last-scattering surface by approximately the 
same amount. This approximation is good at the percent level, 
and can readily be related to the full trispectrum estimator, as dis- 
cussed in more detail in Pearson et al. (2012, Sect. IV); see also 
Okamoto & Hu 2002; Munshi et al. 201 1; Smidt et al. 2010. 

A large-scale power modulation therefore translates directly 
into a large-scale modulation of the small-scale CMB tempera- 
ture: 

T(h) * T g (h)[l + 0(h, r.)] = T g (n)[l + /(«)], (67) 

where T g are the usual small-scale Gaussian CMB temperature 
anisotropies and r* is the radial distance to the last-scattering 
surface. We can quantify the trispectrum as a function of modu- 
lation scale by using the power spectrum of the modulation, 



t nl (L) = — . 

L L 



(68) 



As is conventional, we normalize relative to Cj* , the power spec- 
trum of the primordial curvature perturbation at the location of 
the recombination surface. The field / is directly observable, 
but is not, since the curvature perturbation can only be con- 
strained very indirectly on very large scales. We shall therefore 
give constraints on /, which is directly constrained by Planck, 
but also on tnl for comparison with the inflation literature. Note 
that t nl ~ 500 corresponds to an / = (9(10~ 3 ) modulation. 

A general quadratic estimator methodology for reconstruct- 
ing / was developed in Hanson & Lewis (2009), which we 
broadly follow here. The structure is essentially identical to that 
for lensing reconstruction (Planck Collaboration XVII 2013), 
where here instead of reconstructing a lensing potential (or de- 
flection angle), we are reconstructing a scalar modulation field. 
The quadratic maximum likelihood estimator for the large-scale 
modulation field / (assuming it is small) is given by 



flM - f^LML'M 1 \fL'M' ~ <./z/M'}] 



(69) 



where / is a quadratic function of the filtered data that can be 
calculated quickly in real space: 



Tlm 



Uimi 
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f2>tl2 



(70) 



ax) = &(*)[i + #*)], 



(66) 



Here f' = C l T' is an inverse-variance filtering sky map (which 
accounts for sky cuts and inhomogeneous noise), and CV, is the 
lensed Q. The "mean field" / MF = (/) can be estimated from 
simulations, along with the Fisher normalization T that is given 
by the covariance of / - / MF . The i, j indices are included here, 
since we shall be using different sky maps with independent 
noise to avoid noise biases at the level of modulation field re- 
construction. For low L and high ( mix the reconstruction noise 
is very nearly constant (white, because each small patch of sky 
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gives a nearly-uncorrelated but noisy estimate of the small-scale 
power), and the reconstruction is very local. 

In practice the inverse-variance filtering is imperfect, the 
noise cannot be modelled exactly, and the normalizing Fisher 
matrix Timi'm' evaluated from simulations would be inaccu- 
rate. Instead we focus on / directly, which is approximately an 
inverse-variance weighed reconstruction of the modulation, and 
is manifestly very local in real space (and hence zero in the cut 
part of the sky). Since the reconstruction noise, which also ap- 
proximately determines the normalization, is nearly white (con- 
stant in L), f - / MF in real space has an expectation nearly pro- 
portional to the underlying modulation outside the mask. 

We then define an estimator of the modulation power spec- 
trum 



C 



f 



2L 



M 



fMF\2 
' JLM I 
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where 



iVf = 



2L 



(71) 



(72) 



M 



is a noise bias for zero signal estimated from simulations. The 
normalization Ai is the analytic ideal full-sky normalization 
which is very close to a constant, and ki is a calibration fac- 
tor determined from with-signal simulations. On small scales 
ki oc / sk ', but has some scale dependence: it increases towards 

ki oc / sk ^ at L = at very low L. We shall sometimes plot 

c£ = k~^& L , corresponding to the uncalibrated reconstruction 
of the power modulation, which is very local in real space. 

For each value of the modulation scale L, Eq. (71) defines a 
separate estimator for tnl 



C 



f 



L L 

We can combine estimators from all I by constructing 



(73) 



f N L,l ~ AT 



£ C c * C ov- L } L C f L , 



(74) 



where N = Ez7=l, ^u cov JiPl m ^ cov lu is the covariance of 
from simulations with tnl = 0. On the full-sky the estimators 
from each L would be independent, but the mask introduces sig- 
nificant coupling between the very low multipoles and this form 
of the estimator allows us to account for this. In the full-sky 
uncorrelated approximation, with a nearly scale-invariant pri- 
mordial spectrum and using the whiteness of the reconstruction 
noise, the estimator for tnl simplifies to (Pearson et al. 2012) 



2L+1 C[ 



(75) 



where N = L m j n - (L max + 1) . This result does not require 
many simulations to estimate the covariance accurately for inver- 
sion, and is typically expected to give very similar results with 
an error bar that is less than 10% larger. We calculate both as 
a cross-check, but report results for tnl because it is more ro- 
bust, and in our simulation results actually has slightly lower 
tails (though larger variance). Mean fields and the bias are 
estimated in all cases from 1000 zero-T^L simulations, and the 
mask used retains about 70% of the sky. 



f c 

If there is a nearly scale-invariant signal, so C L oc C L as ex- 
pected in most multi-field inflation models, the contributions fall 
rapidly oc 1/L 3 , as expected when measuring a scale -invariant 
signal that has large white reconstruction noise. The signal is 
therefore on very large scales, with typically half the Fisher sig- 
nal in the dipole modulation and 95% of the signal at L < 4, 
justifying the squeezed approximations used. We use L max = 10 
for the estimators, which includes almost all of the signal-to- 
noise but avoids excessive contamination with the 'blue' spec- 
trum of lensing contributions. However, due to the small number 
of modes involved, the posterior distributions of tnl can have 
quite broad tails corresponding to the finite probability that all 
the largest-scale modulation modes just happen to be near zero. 
To improve constraints on large values it can help to include a 
larger range of L, and we consider L up to L max = 50, which is 
about the limit of where the approximations are valid. 



5. Non-primordial contributions to the CMB 
bispectrum and trispectrum 

In this subsection we present the steps followed to account for 
and remove the main non-primordial contributions to CMB NG. 

5.1. Foreground subtraction 

Foreground emission signals in the microwave bands have a 
strong non-Gaussian signature. Therefore any residual emis- 
sion in the CMB data can give a spurious apparent primordial 
NG detection. In this Section we describe how the foreground 
emission was modelled and treated in the creation of the CMB 
maps used in the present analysis. A comprehensive descrip- 
tion of the sky modelling can be found in Delabrouille et al. 
(2012). The foreground cleaning techniques are described in 
Planck Collaboration XII (2013). 

The pre-launch version of the Planck Sky Model (PSM) is 
based on observations of the emission from our own Galaxy and 
known extra-Galactic sources, largely in the radio and infrared 
bands. The PSM is described in Delabrouille et al. (2012) and 
includes models of CMB (including a dipole), diffuse Galactic 
emissions (synchrotron, free-free, thermal dust, Anomalous 
Microwave Emission and CO molecular lines), emission from 
compact objects (thermal SZ effect, kinetic SZ effect, radio 
sources, infrared sources, correlated far-infrared background and 
ultra-compact Hn regions). The sky model includes total inten- 
sity as well as polarization, which was not used in this paper. 

The PSM has been used to create the sixth round of Full 
Focal Plane (FFP6) simulations, a set of simulations for the 2013 
data release based on detailed models of the sky and instrument 
(e.g., noise properties, beams, satellite pointing and map-making 
process), consisting of both Gaussian and non-Gaussian CMB 
realizations. A detailed description of the FFP6 simulations can 
be found in Planck Collaboration ES (2013). The FFP6 simula- 
tions have been used to test and validate the component separa- 
tion algorithms employed in Planck and to select the ones to be 
applied to the data (Planck Collaboration XII 2013). 

CMB foreground-cleaned maps have been created using 
several independent techniques: explicit parametrization and 
fitting of foregrounds in real space (Commander-Ruler, C-R); 
Spectral Matching of foregrounds implementing Independent 
Component Analysis (SMICA) (Delabrouille et al. 2003; 
Cardoso et al. 2008); Internal Linear Combination (Needlet 
Internal Linear Combination, (NILC) (Delabrouille et al. 
2009); and Internal Template Fitting (Spectral Estimation via 
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Table 1. The bias in the three primordial /nl parameters due to 
the ISW-lensing signal for the four component-separation meth- 
ods. 





SMICA 


NILC 


SEVEM 


C-R 


Local 


7.1 


7.0 


7.1 


6.0 


Equilateral 


0.4 


0.5 


0.4 


1.4 


Orthogonal 


-22 


-21 


-21 


-19 



Table 2. Results for the amplitude of the ISW-lensing bispec- 
trum from the SMICA, NILC, SEVEM, and C-R foreground-cleaned 
maps, for the KSW, binned, and modal (polynomial) estimators; 
error bars are 68% CL . 





SMICA 


NILC 


SEVEM 


C 


-R 


KSW 
Binned . . 
Modal . . . 


0.81 + 0.31 
. . 0.91 ± 0.37 
. . 0.77 ± 0.37 


0.85 ± 0.32 
1.03 ± 0.37 
0.93 ± 0.37 


0.68 ± 0.32 
0.83 + 0.39 
0.60 ± 0.37 


0.75 : 
0.80: 
0.68: 


t0.32 
t0.40 
t0.39 



Expectation Minimization, SEVEM). These and other techniques 
underwent a pre-launch testing phase (Leach et al. 2008). Each 
method provides a Planck CMB foreground-cleaned map with a 
confidence mask, which defines the trusted cleaned region of the 
sky; an estimate of the noise in the output CMB map obtained 
from half-ring difference maps; and an estimate of the beam 
transfer function of the processed map. The resolution reaches 
5 arcminutes. In addition a union of all the confidence masks, 
denotes as U73, is provided. Channels from both the Low 
Frequency Instrument (LFI, Planck Collaboration II 2013) and 
the High Frequency Instrument (HFI, Planck Collaboration VI 
2013) of Planck are used to achieve each of the reconstructed 
CMB templates. The validation of CMB reconstruction through 
component separation is based on the inspection of several 
observables, as explained in detail in Planck Collaboration XII 
(2013): the two-point correlation function and derived cosmo- 
logical parameters; indicators of NG including the /nl results 
presented in the present paper; and cross-correlation with 
known foreground templates. MC simulations varying the CMB 
realizations in the FFP6 simulations were used to establish 
uncertainties on the observables listed above. Based on various 
figures-of-merit, the foreground cleaning techniques performed 
comparably well (Planck Collaboration XII 2013). 

5.2. The Integrated Sachs-Wolfe-lensing bispectrum 

One of the most relevant mechanisms that can generate NG from 
secondary CMB anisotropics is the coupling between weak lens- 
ing and the ISW (Sachs & Wolfe 1967) effect. This is in fact 
the leading contribution to the CMB secondary bispectrum with 
a blackbody frequency dependence (Goldberg & Spergel 1999; 
Verde & Spergel 2002; Giovi et al. 2005). 

Weak lensing of the CMB is caused by gradients in the 
matter gravitational potential that distorts the CMB photon 
geodesies. The ISW on the other hand arise because of time- 
varying gravitational potentials due to the linear and nonlinear 
growth of structure in the evolving Universe. Both the lensing 
and the ISW effect are then related to the matter gravitational 
potential and thus are correlated phenomena. This gives rise to 
a non-vanishing 3-point correlation function. Furthermore, lens- 
ing is related to nonlinear processes which are therefore non- 
Gaussian. A detailed description of the signal, which accounts 



also for the contribution from the early-ISW effect, can be found 
in Lewis (2012). 

The ISW-lensing bispectrum takes the form: 



where P, L, and ISW indicate primordial, lensing and ISW con- 
tributions respectively. This becomes 



nm 1 m 2 m 3 (ISW-L) _ f,m\rmms 1 1SW-L 



(77) 



where Q^l m is the Gaunt integral and M S *7 L is the reduced 
bispectrum given by 
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Here Cj 1 is the lensed CMB power spectrum and C 
is the ISW-lensing cross-power spectrum (Lewis 2012; 
Goldberg & Spergel 1999; Verde & Spergel 2002; Cooray & Hu 
2000) that expresses the statistical expectation of the correlation 
between the lensing and the ISW effect. 

As shown in Hanson et al. (2009b), Mangilli & Verde 
(2009), and Lewis et al. (201 1), the ISW-lensing bispectrum can 
introduce a contamination in the constraints on primordial local 
NG from the CMB bispectrum. Both bispectra are maximal for 
squeezed or nearly squeezed configurations. The bias on a pri- 
mordial /nl (e.g., local) due to the presence of the ISW-lensing 
cross correlation signal is defined as: 
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(79) 



(80) 



where # ISW ~ L and B p refer respectively to the ISW-lensing and 
the primordial bispectrum, and V is defined below Eq. (35). 

The bias in the estimation of the three primordial /nl from 
Planck is given in Table 1 . As one can see, taking into account 
the /nl statistical error bars shown, e.g., in Table 8, the local 
shape is most affected by this bias (at the level of more than 
1 0"iocaiX followed by the orthogonal shape (at the level of about 
0.5 cr ort ho), while the equilateral shape is hardly affected. In this 
paper we have taken into account the bias reported in Table 1 by 
subtracting it from the measured /nl- 7 

The results for the amplitude of the ISW-lensing bispectrum 
from the different foreground-cleaned maps are given in Table 2. 
It should be noted that the binned and modal estimators are 
less correlated to the exact template for the ISW-lensing shape 
than they are for the primordial shapes, hence their larger er- 
ror bars compared to KSW (which uses the exact template by 
construction (Mangilli et al. 2013)). The conclusion is that we 
detect the ISW-lensing bispectrum at a value consistent with the 
fiducial value of 1, at a significance level of 2.6 cr (taking the 
SMICA-KSW value as reference). For details about comparisons 
between different estimators and analysis of the data regarding 
primordial shapes we refer the reader to Sects. 6 and Sect. 7. 



See Kim et al. (2013) for other debiasing techniques. 
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Fig. 2. The binned skew-Q statistics from the SMICA map for (a) 
ISW-lensing and (b) Poisson point sources. Theoretical curves 
are not fitted to the data shown, but are plotted with the ampli- 
tude (the only free parameter) determined from the KSW tech- 
nique. The Poisson point-source foreground is clearly detected, 
and the ISW-lensing skew-spectrum is evident for I < 1750, 
with a suggestion of another source of NG at high {, b ps is the 
Poisson point-source amplitude in dimensionless units of 10~ 29 , 
and /^l W L i s me ISW-lensing amplitude in units of that expected 
from the Planck best-fit cosmology. Note that error bars are from 
data-averaging, and as a consequence are underestimates. 



We show for the SMICA map in the top figure of Fig. 2 
the measured skew-CV spectrum (see Sect. 3.2.2) for opti- 
mal detection of the ISW-lensing bispectrum, along with the 
best-fitting estimates of /nl from the KSW method for dif- 
ferent values of I. It should be noted that the skew-Q spec- 
trum is not a fit to the KSW data points; its shape is fully 
fixed by the template under consideration, with only the over- 
all amplitude as a free parameter. Hence the agreement be- 
tween the curve and the points in the regime up to I ~ 
1750 is good evidence that KSW is really detecting the ISW- 
lensing effect and not some other source of NG (although there 
might be some evidence of an additional NG contribution at 
{ > 1750; note that point sources, at the level determined 
by their own skew-spectrum, do not contribute significantly 
to the ISW-lensing statistic). See Planck Collaboration XVII 
(2013), Planck Collaboration XIX (2013) for further informa- 
tion about the detection by Planck of the ISW-lensing signal. 



5.3. Point-sources bispectrum 

Extra-Galactic point sources at Planck frequencies are divided 
into two broad categories: radio sources with synchrotron and/or 
free-free emission; and infrared galaxies with thermal emission 
from dust heated by young stars. Radio sources are dominant at 
central CMB frequencies up to 143 GHz, and can be considered 
unclustered (Toffolatti et al. 1998; Gonzalez-Nuevo et al. 2005). 
Hence their bispectrum is constant and is related to their number 



Table 3. Results for the amplitude of the point source (Poisson) 
bispectrum (in dimensionless units of 10~ 29 ) from the SMICA, 
NILC, SEVEM, and C-R foreground-cleaned maps, for the KSW, 
binned, and modal (polynomial) estimators; error bars are 68% 
CL. Note that the KSW and binned estimators use ^ max = 2500, 
while the modal estimator has ^ max = 2000. 





SMICA 


NILC 


SEVEM 


C- 


-R 


KSW 


7.7 + 1.5 


9.2 ± 1.7 


7.6 ± 1.7 


1.1 


± 5.1 


Binned . . 


. 7.7 ±1.6 


8.2 ± 1.6 


7.5 ± 1.7 


0.9 


+ 4.8 


Modal . . . 


. 10 + 3 


11 ±3 


10 ±3 


0.5 


± 6 



counts as 



J* 

Jo 



5 dn 
dS' 



'rdS, 



(81) 



with S the flux density, dn/dS the number counts per steradian, 
S c the flux cut and k v the conversion factor from flux to relative 
temperature elevation, depending on the frequency and instru- 
mental bandpass. 

Infrared galaxies become important at higher frequencies, 
217GHz and above, and are highly clustered in dark matter ha- 
los, which enhances their bispectrum on large angular scales 
(Lacasa et al. 2012; Curto et al. 2013). However, in the Planck 
context it was shown by Lacasa & Aghanim (2012) that the IR 
bispectrum is more than 90% correlated with the Poissonian 
template of the radio sources. So a joint estimation of /nl with a 
Poissonian bispectrum template will essentially account for the 
IR signal, and provide quasi-identical values compared to an 
analysis accounting for the IR bispectrum template. Indeed, in 
our final optimal bispectrum constraints for primordial shapes, 
we will account for the potential contamination from point 
sources by jointly fitting primordial and Poisson templates to the 
data. 

Our final measured point-source bispectrum amplitudes 
from the data are reported in Table 3. The amplitude is expressed 
in dimensionless units, i.e., it has been divided by the appropri- 
ate power of the monopole temperature To, and has been mul- 
tiplied by 10 29 . As shown in Sect. 8.1, the Poisson template is 
the only one that still evolves significantly between I = 2000 
and i = 2500. This explains the differences between the values 
of the KSW and binned (that use ( max - 2500) and the modal 
(that uses £ max = 2000) estimators. It has been shown that for 
the same value of i max all three estimators agree very well. 

We finally conclude from Table 3 that we detect the point- 
source bispectrum with high significance in the SMICA, NILC, 
and SEVEM cleaned maps, while it is absent from the C-R cleaned 
map. The measured skew-Q spectrum of the SMICA map in the 
bottom figure of Fig. 2 gives further evidence that the NG from 
foreground point sources is convincingly detected. The only de- 
gree of freedom in this plot is the amplitude, which is not set by 
a direct fit to the skew-Q , but rather is estimated by KSW. As a 
result, the good agreement with the shape of this skew-Q spec- 
trum is powerful evidence that there is NG from point sources. 
However, this still turns out to be a negligible contaminant for 
primordial /nl studies, due to the very low correlation between 
the Poisson bispectrum and the primordial shapes. 

5.4. Non-primordial contributions to the trispectrum 

The main non-instrumental source of non-primordial sig- 
nal is the kinematic modulation dipole due to the pecu- 
liar velocity of the earth, v, whose magnitude is 0(1 0~ 3 ) 
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(Challinor & van Leeuwen 2002; Kosowsky & Kahniashvili 
2011; Amendola et al. 2011). If data are used to constrain 
tnl using the dipole modulation (which shrinks the Fisher 
error by a factor of two relative to starting at L — 2), the 
dipole-induced signal must be subtracted, since its modulation 
reconstruction has signal-to-noise larger than unity at Planck 
resolution. Confirmation that this signal is detected with the 
expected magnitude and direction is a good test of our method- 
ology. The dipole signal seen by Planck is studied in detail 
in Planck Collaboration XXVII (2013), so we only summarize 
the key points here. 

The local Doppler effect modulates the observed CMB tem- 
perature Tq[1 + AT(h)] by 1 + h ■ v at leading order, so that 
T(n) = 7o[l + Ar(«)](l + h ■ u). The spectrum in each direction 
remains a blackbody, but the relative response in the intensity 
I v (v, h) at the observed Planck frequencies is however frequency 
dependent. The effective thermodynamic fractional temperature 
anisotropy A© at each frequency for zero peculiar velocity is 
defined by 



I v (v, h) = / y (v) 



1 + 



dln/ v 



dlnT 



A©(«) 



(82) 



With peculiar velocity the temperature T depends on the second 
order term ATh ■ v, so expanding the Planck function to second 
order then gives a change in the effective small-scale temper- 
ature anisotropy from both first and second order terms in the 
expansion of I v : 



A©(«) 



1 + n ■ v + T— — — — n ■ v 



A@(n) 



dIJdT 

(1 + [xcoth(x/2) - 1] h ■ v)A®(h), 



(83) 



where x = hv/k B To (and we neglect small second-order non- 
modulation terms). Thus the anisotropies in the Planck maps 
have a dipolar modulation given by 1 + b v h ■ v, where for the fre- 
quency bands we use ~ 2, £>2i7 ~ 3, and/? = \v\ = 1.23xl0 -3 
in the direction of CMB dipole. In addition our peculiar velocity 
induces kinetic aberration, which looks at leading order exactly 
like a dipole lensing convergence and only projects weakly into 
the power anisotropy estimator. For constraining tnl both of the 
expected kinematic signals can be included in the simulations, 
and hence subtracted in the mean field of the modulation recon- 
struction. 

Secondary effects are dominated by the significant and very 
blue lensing signal. However unlike for the /nl bispectmm lens- 
ing only overlaps with tnl at a small fraction of the error bar as 
long as only low modulation multipoles L < 10 are used, where 
the tnl signal peaks (Pearson et al. 2012). We include lensing in 
the simulations, so lensing is straightforwardly accounted for in 
our analysis by its inclusion in the A^ 0) noise bias (Eq. (72)) and 
mean field. 

A variety of instrumental effects can also give a spurious 
modulation signal if not modelled accurately. In particular the 
mean field due to anisotropic noise is very large (Hanson et al. 
2009a). On the ultra-large scales of interest for tnl, our un- 
derstanding of the noise is not adequate to calculate accurately 
and subtract this large signal. Instead, as for the power spec- 
trum estimation, we use cross-map estimators that have no noise 
mean field on average. Both the noise and most other instru- 
mental effects such as gain variations are expected to produce 
a signal with approximate symmetry about the ecliptic plane. 
Our modulation reconstruction methodology is especially use- 
ful here, since we can easily inspect the orientation of any sig- 
nal found; for example a naive treatment of the noise not using 



cross-maps would give a large apparent quadrupolar modulation 
signal aligned with the ecliptic, corresponding to percent-level 
misestimation of the noise mean field from inaccurate noise sim- 
ulation. 

Beam asymmetries are included in the simulation, as de- 
scribed in Planck Collaboration XVII (2013), but their effect 
is very small, since the modulation we are reconstructing is 
isotropic. 

Since we are reconstructing a modulation of small-scale 
power, the estimator is totally insensitive to smooth large- 
scale foregrounds. However large-scale variation in small- 
scale foreground power can mimic a trispectrum modulation. 
We project out 857 GHz as a dust template in our inverse- 
variance filtering procedure, as described in Appendix A of 
Planck Collaboration XVII (2013), but do not include any other 
foreground model in the trispectrum analysis. Any unmodelled 
foreground power variation would increase the tnl signal, so our 
modelling is sufficient to place a robust upper limit. 

6. Validation Tests 

The /nl results quoted in this paper have all been cross- 
validated using multiple bispectrum-based estimators from dif- 
ferent groups. Having multiple estimators was extremely useful 
for the entire analysis, for two main reasons. First, it allowed 
great improvement in the robustness of the final results. In the 
early stages of the work the comparison between different in- 
dependent techniques helped to resolve bugs and other techni- 
cal issues in the various computer codes, while during the later 
stages it was very useful to understand the data and find the opti- 
mal way of extracting information about the various bispectmm 
templates. Secondly, besides these cross-checking purposes, dif- 
ferent estimators provide also interesting complementary infor- 
mation, going beyond simple /nl estimation. For example, the 
binned and modal estimators provide a reconstruction of the 
full bispectmm of the data (smoothed in different domains), the 
skew-Q estimator allows monitoring of the contribution to /nl 
from different sources of NG, the wavelets reconstruction allows 
/nl directionality tests, and so on. 

In this Section we are concerned with the first point above, 
that is, the use of multiple bispectrum-based pipelines as a way 
to improve the robustness of the results. For this purpose, a large 
amount of work was dedicated to the development and analy- 
sis of various test maps, in order to validate the estimators. This 
means not only checking that the various estimators recover the 
input /nl within the expected errors, but also that the results 
agree on a map-by-map basis. 

The Section is split into two parts. Sect. 6. 1 shows results on 
a set of initially full-sky, noiseless, Gaussian CMB simulations, 
to which we add, in several steps, realistic complications, in- 
cluding primordial NG, anisotropic coloured noise, and a mask, 
showing the impact on the results at each step. In Sect. 6.2 we 
show our results on a set of simulations that mimic the real data 
as closely as possible (except for the presence of foreground 
residuals, which will be studied in Sect. 8.4): no primordial NG, 
but NG due to the ISW-lensing effect; simulated instrumental 
effects and realistic noise; and simulations passed through the 
component separation pipelines. In fact these are the simulations 
that are used to determine the error bars for the final Planck re- 
sults. 

We present here only a small subset of the large number of 
validation tests that were performed. For example, we also had a 
number of "blind /nl challenges", in which the different groups 
received a simulated data set with an unknown value of input /nl 
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for a given shape and they had to report their estimated values. In 
addition different noise models were tested (white vs. coloured 
and isotropic vs. anisotropic), leading to the conclusion that it 
is important to make the noise in the estimator calibration as 
realistic as possible (coloured and anisotropic). We also tested 
different Galactic and point source masks, with and without in- 
painting, concluding that it is best to fill in both the point sources 
and the Galactic mask, using a sufficient number of iterations in 
our diffusive procedure to entirely fill in the point source gaps, 
while at the same time only effectively apodizing the Galactic 
mask (no small-scale structure in its interior). There were also 
various tests on realistic simulations of Planck data, including 
detailed modelling of the Planck satellite, and the sky signals 
(Gaussian or non-Gaussian CMB and all foregrounds, provided 
by the PSM). These simulations were tested both before and af- 
ter they had passed through the component separation pipelines. 
In all comparison tests the results were consistent with input /nl 
values and differences between estimators were consistent with 
theoretical expectations. 

6.1. Validation of estimators in the presence of primordial 
non-Gaussianity 

The aim of the first set of validation tests is threefold. First, we 
want to study the level of agreement from different estimators 
in ideal conditions (i.e., full-sky noiseless data). The expected 
scatter between measurements is, in this case, entirely due to 
the slightly imperfect correlation between weights of estimators 
that adopt different schemes to approximate the primordial shape 
templates. For this case the scatter can be computed analytically 
(see Appendix A for details). We can then verify that our results 
in ideal conditions match theoretical expectations. This is done 
in Sect. 6.1.1. Second, we want to make sure that the estima- 
tors are unbiased and correctly recover /nl in input for local, 
equilateral, and orthogonal shapes. This is done in Sects. 6.1.2 
and 6.1.3, where a superposition of local, equilateral and orthog- 
onal bispectra is included in the simulations and the three /nl 
values are estimated both independently and jointly. Finally we 
want to understand how much the agreement between pipelines 
in ideal conditions is degraded when we include a realistic cor- 
related noise component and a sky cut, thus requiring the intro- 
duction of a linear term in the estimators in order to account 
for off-diagonal covariance terms introduced by the breaking 
of rotational invariance. Since we want to study the impact of 
adding noise and masking separately, we will first work on a set 
of full-sky maps with noise in Sect. 6. 1 .2, and then add a mask 
in Sect. 6.1.3. 

The tests that we are going to show were applied to the KSW, 
binned and modal estimators. These are three optimal bispec- 
trum pipelines used to analyse Planck data in Sect. 7. Our goal 
for this set of tests is not so much to attain the tightest possi- 
ble agreement between methods, as it is to address the points 
summarized in the above paragraph. For this reason the estima- 
tor implementations used in this specific Section were slightly 
less accurate but faster to compute than those adopted for the fi- 
nal data analysis of Sect. 7. The primary difference with respect 
to the main analysis is that a smaller number of simulations was 
used to calibrate the linear term (80-100 in these tests, as against 
200 or more for the full analysis). For the modal estimator we 
also use a faster expansion with a smaller number of modes: 300 
here versus 600 in the high accuracy version of the pipeline 8 



8 While most of the modal results in this paper come from the most 
accurate 600 modes pipeline, a few computationally intensive data vali- 



used in Sect. 7. Even with many fewer modes, the modal estima- 
tor is still quite accurate: the correlation coefficient for the modal 
expansion of the local template is 0.95, while for the equilateral 
and orthogonal shapes it is 0.98. 

6.1.1. Ideal Gaussian simulations 

As a basis for the other tests we start with the ideal case, a set 
of 96 simulations of a full-sky Gaussian CMB, with a Gaussian 
beam with FWHM 5 arcmin and without any noise, cut off at 
f m „ = 2000 in our analyses. The independent Fisher matrix er- 
ror bars in that case are 4.2 for local NG, 56 for equilateral, and 
28 for orthogonal. 

Note that this test does not make sense for all estimators, and 
hence results are not included for all of them. For example, for 
the binned estimator the optimal binning depends on the noise. 
While this dependence is not very strong, the difference between 
no noise and Planck noise is sufficiently large that a completely 
different binning would have to be used just for this test, going 
against the purpose of this Section to validate the estimators as 
used for the data analysis. 9 

The purpose here is mostly aimed at checking consistency 
with the following formula (derived in Appendix A) for the ex- 
pected scatter (standard deviation) between / N l results of the 
same map from an exact and an approximate estimator: 

Vl -r 2 

£T 5/NL =Ath . (84) 

Here A t h is the standard deviation of the exact estimator and r 
is the correlation coefficient that gives the correlation of the ap- 
proximate bispectrum template with the exact one, defined as 

E f l f 2 f 3 [ \ [ 2 [ 1 
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where the label "th" denotes the initial bispectrum shape to fit to 
the data, and "exp" is the approximate expanded one. Note that 
this formula has been obtained under the simplifying assump- 
tions of Gaussianity, full-sky coverage and homogeneous noise. 
For applications dealing with more realistic cases we might ex- 
pect the scatter to become larger, while remaining qualitatively 
consistent. 

The results averaged over the whole set of maps are given in 
Table 4 for the KSW and modal estimators individually, as well 
as for their difference. The plane wave modal expansion imple- 
mented here achieves about 98% correlation with the separable 
shapes used by KSW. According to the formula above we then 
expect a standard deviation of map-by-map differences of order 
0.2 A fa for a given shape, where A fa is the corresponding /nl 
error bar. Looking at the left-hand side of Table 4, we see that 
the error bars are 4 for local NG, 50 for equilateral, and 30 for 
orthogonal. So we predict a standard deviation of map-by-map 

dation tests of Sect. 8 also use the fast 300 modes version; therefore the 
results in this Section also provide a direct validation of the fast modal 
pipeline. 

9 While the binning with 48 bins and C max = 2000 used in the vali- 
dation tests of Sects. 6.1.2 and 6.1.3 is also slightly different from the 
binning used for the data analysis with 5 1 bins and f raax = 2500, these 
differences are very small and the binnings have very similar correlation 
coefficients of 0.99 or more for local and equilateral shapes, and 0.95 
for orthogonal. 
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Table 4. Results for /nl for the set of ideal Gaussian simulations 
described in Sect. 6.1.1 for the KSW and modal estimators and 
for their difference, assuming all shapes to be independent. 



Independent 


KSW 


Modal 


Modal - 


KSW 


Local 


-0.5 ±4.1 


-0.5 ±4.1 . 


-0.0 ± 


0.6 


Equilateral 


2.2 ± 48 


1.3 ±48 . 


-0.9 ± 


8.9 


Orthogonal 


-1.1 ± 29 


-1.0 ±30 . 


0.1 ± 


6.5 



differences of 0.8, 10 and 6 for local, equilateral, and orthogo- 
nal NG, respectively. As one can see from the "Modal-KSW" 
column, the measurements are in excellent agreement with the 
theoretical expectation. 

6.1 .2. Non-Gaussian simulations with realistic noise 

A set of 96 full-sky non-Gaussian CMB simulations was created 
according to the process described by Fergusson et al. (2010a), 
with local f^ al = 12, equilateral f^} = 35, and orthogonal 
/°" ho = -22. The effect of a 5 arcmin beam was added, as well 
as realistic coloured and anisotropic noise according to the speci- 
fications of the SMICA cleaned map. The independent Fisher ma- 
trix error bars in that case are 5.3 for local, 63 for equilateral, 
and 33 for orthogonal NG, while the joint ones are respectively 
6.0, 64, and 37. 

The results averaged over the whole set are given in Table 5 
for the various estimators individually, as well as for the differ- 
ences with respect to KSW. Compared to the previous case we 
now deviate from the exact theoretical expectation for two rea- 
sons: we include a realistic correlated noise component; and we 
have NG in the maps. The presence of NG in the input maps 
will lower the agreement between estimators with respect to the 
Gaussian case if the correlation between weights is not exactly 
100%. This is even more true in this specific case, where NG 
of three different kinds is present in the input maps and also 
cross-correlation terms between different expanded shapes are 
involved (and propagated over in the joint analysis). Moreover, 
when noise is included the specific modal expansion used for 
this test is 95% correlated to the separable KSW local shape (so 
there is a 3% reduction of the correlation compared to the ideal 
case for the modal local shape); we thus expect a further degra- 
dation of the level of agreement for this specific case. Finally, in 
order to correct for noise effects, a linear term has to be added 
to the estimators. Since the linear term is obtained by MC aver- 
aging over just 80 or 96 simulations in this test (depending on 
the estimator), MC errors are also adding to the measured differ- 
ences. Of course the MC error can be reduced by increasing the 
number of simulations in the linear term sample. We do this for 
the analysis of the real data and in Sect. 6.2, but it was computa- 
tionally too expensive for this set of preliminary validation tests, 
so we decided here to just account for it in the final interpretation 
of the results. 

As a consequence of the above, we can no longer expect the 
map-by-map /nl differences to follow perfectly the theoretical 
expectation, obtained in the previous Section in idealized con- 
ditions (full-sky, no noise, and Gaussianity). With these caveats 
in mind, the agreement between different pipelines remains very 
good, being about 0.3 cr in most cases and about 0.5 cr for the 
modal-KSW difference in the local case, which can be easily ex- 
plained by the fact that this is the set of weights with the lowest 



correlation (95%, as mentioned above). All estimators are unbi- 
ased and recover the correct input values. 

6.1 .3. Impact of the mask 

To the simulations of Sect. 6.1.2 we now apply the Planck 
union mask - denoted U73 - masking both the Galaxy and 
the brightest point sources and leaving 73% of the sky un- 
masked (Planck Collaboration XII 2013). This is the same mask 
used to analyse Planck data in Sect. 7. The independent Fisher 
matrix error bars in that case (taking into account the / s w cor- 
rection) are 6.2 for local NG, 74 for equilateral, and 39 for or- 
thogonal, while the joint ones are respectively 7.1, 76, and 44. 

All masked areas of the sky (both Galactic and point sources) 
are filled in with a simple iterative method. In this simple in- 
painting method each pixel in the mask is filled with the average 
of all eight surrounding pixels, and this is repeated 2000 times 
over all masked pixels. The filling-in helps to avoid propagating 
the effect of a sharp edge and the lack of large-scale power inside 
the mask to the unmasked regions during harmonic transforms. 
This inpainting method is the one that was used to produce all 
NG results in this paper for methods that need it (KSW, binned 
and modal). 

The results averaged over the whole set of simulations are 
given in Table 6 for the various estimators individually, as well 
as for the differences with respect to KSW. The map-by-map 
results are shown in Fig. 3. 

This is the most realistic case we consider in this set of tests. 
Besides noise, we also include a sky cut and our usual mask in- 
painting procedure. All the caveats mentioned for the previous 
case are still valid, and possibly emphasized by the inclusion of 
mask and inpainting. In the light of this, the agreement is still 
very good, worsening a bit with respect to the "full-sky + noise" 
case only for the local measurement, where the mask is indeed 
expected to have the biggest impact. In the joint analysis all esti- 
mators recover the correct input values for the local and orthog- 
onal cases, but all estimators find a value for equilateral NG that 
is somewhat too low. It is unclear whether this is an effect of 
masking and inpainting on the equilateral measurement or just a 
statistical fluctuation for this set of simulations. In any case, this 
potential bias is small compared to the statistical uncertainty, so 
that it would not have a significant impact on the final results. 

To summarize the results of this Sect. 6.1, we performed an 
extensive set of validation tests between different /nl estimators 
using strongly, but not perfectly, correlated primordial NG tem- 
plates in their weights. The test consisted in comparing the / NL 
measured by the different estimators for different sets of simula- 
tions, on a map-by-map basis. We started from ideal conditions: 
full-sky Gaussian noiseless maps. In this case we computed a 
theoretical formula providing the expected standard deviation of 
the /nl differences, as a function of the correlations between the 
input NG templates in the different estimators. Our results match 
this formula very well. In the other two simulation sets we added 
realistic features (noise, mask and inpainting) and we included 
a linear combination of local, equilateral and orthogonal NG. 
First of all we verified that all the pipelines correctly recover 
the three /nl input values, hence they are unbiased. Moreover, 
we observed that adding such features produces an expected 
slight degradation of the level of agreement between different 
pipelines, that nevertheless remains very good: about 0.3-0.4 cr 
for equilateral and orthogonal NG, and about 0.5-0.6 cr for local 
NG, which is the shape most affected by mask and noise con- 
tamination. 
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Table 5. Results from the different estimators for /nl for the set of full-sky non-Gaussian simulations described in Sect. 6.1.2. Both 
the results for the estimators individually and for the differences with KSW are given. 





KSW 


Binned 


Modal 


Binned - 


- KSW 


Modal - 


- KSW 


Independent 
















Local 


13.8 + 5.2 


14.1 ± 5.2 


14.1 + 5.3 


0.3 + 


2.1 


0.4 + 


2.6 


Equilateral 


63 + 57 


62 ±58 


64 + 57 


-0.9 + 


20.1 


1.0 + 


18.1 


Orthogonal 


-52 + 37 


-58 ± 40 


-54 + 37 


-6.0 + 


12.6 


-2.2 + 


12.3 


Joint 
















Local 


11.7 + 6.2 


12.0 + 6.6 


12.0 + 6.4 


0.2 + 


2.7 


0.2 ± 


3.2 


Equilateral 


31 + 59 


29 + 61 


31 + 59 


-1.8 + 


21.1 


-0.2 + 


18.5 


Orthogonal 


-20 + 43 


-22 ± 47 


-21 +42 


-2.1 + 


15.6 


-0.6 + 


14.8 



Table 6. Results from the different estimators for /nl for the set of masked non-Gaussian simulations described in Sect. 6.1.3. Both 
the results for the estimators individually and for the differences with KSW are given. 





KSW 


Binned 


Modal 


Binned - 


KSW 


Modal - 


- KSW 


Independent 
















Local 


13.5 + 7.1 


13.1 + 6.5 


14.0 + 6.8 


-0.3 + 


3.5 


0.5 + 


4.6 


Equilateral 


55 + 64 


50 + 59 


58 + 63 


-4.4 + 


24.1 


3.3 + 


20.2 


Orthogonal 


-50 + 45 


-53 ± 46 


-52 + 45 


-3.5 + 


16.4 


-1.9 + 


15.2 


Joint 
















Local 


11.7 + 8.3 


11.4 + 7.9 


12.2+8.4 


-0.3 + 


4.3 


0.4 + 


5.7 


Equilateral 


23 + 66 


19 + 59 


24 + 64 


-3.8 + 


27.7 


1.7 ± 


24.8 


Orthogonal 


-18 + 51 


-20 ± 54 


-18 + 55 ... 


-1.3 + 


19.9 


0.3 ± 


20.4 



Table 7. Results from the different estimators for /nl for 99 maps from a set of realistic lensed simulations passed through the 
SMICA pipeline, described in Sect. 6.2. Both the results for the estimators individually and for the differences with KSW are given. 



Independent 


KSW 


Binned 


Modal 


Wavelet 


Binned - 


KSW 


Modal - 


- KSW 


Wavelet 


- KSW 


Local 


7.6 + 6.0 


6.8 + 5.8 


7.7 + 5.9 


8.1 + 8.4 


-0.8 + 


1.2 


0.1 


± 1.4 


0.5 ± 


6.4 


Equilateral 


4 + 76 


-1 ± 72 


2 + 76 


-3 + 76 


-5 + 


20 


-2 


± 13 


-7 ± 


91 


Orthogonal 


-21 + 42 


-20 ± 41 


-21+42 


-15 + 53 ... 


1.6 + 


11 


-0.1 


± 8 


6.4 + 


48 



6.2. Validation of estimators on realistic Planck simulations 

In the tests of the previous Subsection we checked the bias of the 
estimators and studied their level of agreement, given the corre- 
lation between their weights, in the presence of noise and a sky 
cut. To speed up the computation while still retaining enough 
accuracy for the purposes of that analysis, we used a relatively 
small number of maps for linear term calibrations (80-100) and 
used a smaller number of modes than usual in the modal esti- 
mator. In the present Subsection we instead try to simulate as 
accurately as possible real data analysis conditions. Our goal is 
to obtain an accurate MC-based expectation of the scatter be- 
tween different /nl measurements when the pipelines are run on 
actual Planck maps. 

To this aim we use FFP6 simulation maps described in 
Planck Collaboration ES (2013). The original FFP6 maps were 
lensed using the Lenspix algorithm, and processed through the 
SMICA component separation pipeline. They were then multi- 
plied by the Galactic and point source mask U73 as in the ac- 
tual /nl analysis, and inpainted as usual. Since our final re- 
sults show full consistency with Gaussianity for local, equilateral 
and orthogonal shapes, we do not include any primordial / N l in 
these maps. We note that although the simulations were passed 
through SMICA in order to provide a realistic filtering of the data, 



they did not include any foreground components. The impact of 
foreground residuals will be studied separately in Sect. 8.4. 

The configuration of all bispectrum pipelines was the same 
as used for the final data analysis, which implies a correlation 
of 99% or better between the weights of the KSW, binned and 
modal estimators. Linear terms were calibrated using 200 sim- 
ulations, after verifying that this number allows accurate con- 
vergence for all the shapes. For this test we also included the 
wavelet bispectrum pipeline described in Sect. 3. Although this 
last estimator turns out to be about 30% suboptimal and, in its 
current implementation, less correlated with the primordial tem- 
plates than the other estimators, it does provide an additional 
interesting cross-check of our results by introducing another de- 
composition basis. We thus used it to analyse SMICA data in 
Sect. 7, while the other three pipelines were used on all maps. 

A comparison of the measured /nl map-by-map for all 
shapes and estimators is shown in Fig. 4. As an overall figure of 
merit of the level of agreement achieved by different pipelines 
we take as usual the standard deviation of the map-by-map /nl 
differences, o~s f . Table 7 shows that the final agreement be- 
tween the three optimal pipelines (KSW, binned, and modal) is 
close to saturating the ideal bound in Eq. (84) determined by the 
imperfect correlation of the weights, i.e., it varies from about 
once to twice cr Sf ^ ^ 0.15 A/ NL for an r = 0.99 correlation. This 
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Fig. 3. Map-by-map comparison of the results from the differ- 
ent estimators for local (top), equilateral (centre), and orthogo- 
nal (bottom) /nl for the set of masked non-Gaussian simulations 
described in Sect. 6.1.3, assuming the shapes to be independent. 
The horizontal solid line is the average value of all maps for 
KSW, and the dashed and dotted horizontal lines correspond to 
lcr and 2cr deviations, respectively. 
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Fig. 4. Map-by-map comparison of the results from the differ- 
ent estimators for local (top), equilateral (centre), and orthogonal 
(bottom) /nl for 99 maps from a set of realistic lensed simula- 
tions passed through the SMICA pipeline, described in Sect. 6.2, 
assuming the shapes to be independent. The horizontal solid line 
is the average value of the maps for KSW, and the dashed and 
dotted horizontal lines correspond to lcr and 2<x deviations, re- 
spectively. 



is very consistent with the level of agreement that we find be- 
tween estimators for the final results from the data, providing a 
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good indication that no spurious NG features are present in the 
actual data set when compared to our simulations. It should be 
noted that we found a similarly good level of agreement between 
estimators for the non-primordial shapes of point sources and 
ISW-lensing, although we chose not to present those results here 
in order to focus on the primordial shapes. Finally, regarding the 
wavelet pipeline, the lower weight correlation and suboptimal 
error bars produce an expected larger scatter when compared to 
the other estimators. Nonetheless, the level of agreement is still 
of order 1 cr, which is quite acceptable for consistency checks of 
the optimal results. Again, this MC expectation agrees with what 
we see in our results on the real data. 

7. Results 

For our analysis of Planck data we considered foreground- 
cleaned maps obtained with the four component separation 
methods SMICA, NILC, SEVEM, and C-R. For each map, / NL 
amplitudes for the local, equilateral, and orthogonal primordial 
shapes have been measured using three (four for SMICA) bispec- 
trum estimators described in Sect. 3. The results can be found 
in Sect. 7.1. These estimators, as explained earlier, basically use 
an expansion of the theoretical bispectrum templates in different 
domains, and truncate the expansion when a high level of corre- 
lation with the primordial templates is achieved. These accurate 
decompositions, which are highly correlated with each other, are 
then matched to the data in order to extract /nl- The different 
expansions are all different implementations of the maximum- 
likelihood estimator given in Eq. (32). So the final estimates are 
all expected to be optimal, and measure /nl from nearly identi- 
cal fitting templates. As discussed and tested in detail on simu- 
lations in Sect. 6, central /nl values from different methods are 
expected to be consistent with each other within about 0.3cry NL . 
It is then clear that comparing outputs from both different esti- 
mators and different component separation methods, as we do, 
allows for stringent internal consistency checks and improved 
robustness of the final / NL results. 

In addition, the binned and modal techniques produce shape- 
independent full bispectrum reconstructions in their own differ- 
ent domains. These reconstructions, discussed in Sect. 7.2, com- 
plement the standard /nl measurements in an important way, 
since they allow detection of possible NG features in the three- 
point function of the data that do not correlate significantly with 
the standard primordial shapes. This advantage is shared by the 
skew-Q method, also applied to the data. A detection of such 
features would either produce a warning that some residual spu- 
rious NG effects are still present in the data or provide an in- 
teresting hint of "non-standard" primordial NG that is not cap- 
tured by the local, equilateral and orthogonal shapes. Additional 
constraints for a broad range of specific models are provided 
in Sect. 7.3 (see also Sect. 2.3). In Sect. 7.4 we present the 
constraints on local NG obtained with Minkowski Functionals. 
Finally, in Sect. 7.5 we present our CMB trispectrum results. 

7.1. Constraints on local, equilateral and orthogonal / N l 

Our goal here is to investigate the standard separable local, equi- 
lateral and orthogonal templates used e.g., in previous WMAP 
analyses (see e.g., Bennett et al. 2012). When using the modal, 
binned, or wavelet estimator, these theoretical templates are ex- 
panded approximately (albeit very accurately) using the relevant 
basis functions or bins. On the other hand, the KSW estimator by 
construction works with the exact templates and, for this reason, 
it is chosen as the baseline to provide the final /nl results for 



Table 8. Results for the /nl parameters of the primordial local, 
equilateral, and orthogonal shapes, determined by the KSW es- 
timator from the SMICA foreground-cleaned map. Both indepen- 
dent single-shape results and results marginalized over the point 
source bispectrum and with the ISW-lensing bias subtracted are 
reported; error bars are 68% CL . 





Independent 


ISW-lensing subtracted 




KSW 


KSW 


SMICA 






Local 


9.8 + 5.8 


2.7 ± 5.8 


Equilateral .... 


-37 + 75 


-42 ± 75 


Orthogonal . . . . 


-46 + 39 


-25 ± 39 



the standard shapes (local, equilateral, orthogonal), see Table 8. 
However, both the binned and modal estimators achieve opti- 
mal performance and an extremely high correlation for the stan- 
dard templates (~ 99%), so they are statistically equivalent to 
KSW, as demonstrated in the previous section. This means that 
we can achieve a remarkable level of cross-validation for our 
Planck NG results. We will be able to present consistent con- 
straints for the local, equilateral and orthogonal models for all 
four Planck foreground-cleaned maps, using three independent 
optimal estimators (refer to Table 9). Regarding component sep- 
aration methods, we adopt the SMICA map as the default for the 
final KSW results given its preferred status among foreground- 
separation techniques in Planck Collaboration XII (2013). The 
other component separation maps will be used for important 
cross-validation of our results and to evaluate potential sensi- 
tivity to foreground residuals. 

All the results presented in this Section were obtained using 
the union mask U73, which leaves 73% of the sky unmasked. 
The mask is the union of the confidence masks of the four differ- 
ent component separation methods, where each confidence mask 
defines the region where the corresponding CMB cleaning is 
trusted (see Planck Collaboration XII 2013). As will be shown in 
Sect. 8.2, results are robust to changes that make the mask larger, 
but choosing a significantly smaller mask would leave some NG 
foreground contamination. For the linear term CMB and noise 
calibration, and error bar determination, we used sets of realistic 
FFP6 maps that include all steps of data processing, and have 
realistic noise and beam properties (Planck Collaboration ES 
2013). The simulations were also lensed using the Lenspix al- 
gorithm and filtered through the component separation pipelines. 

In Table 8 we show results for the combination of the KSW 
estimator and the SMICA map, at a resolution of ( max - 2500. 
We present both "independent" single-shape results and "ISW- 
lensing subtracted" ones. The former are obtained by directly 
fitting primordial templates to the data. For the latter, two ad- 
ditional operations have been performed. In the first place, as 
the name indicates, they have been corrected by subtracting 
the bias due to the correlation of the primordial bispectra to 
the late-time ISW-lensing contribution (Mangilli & Verde 2009; 
Junk & Komatsu 2012; Hanson et al. 2009b, see Sect. 5.2). In 
addition, a joint fit of the primordial shape with the (Poissonian) 
point source bispectrum amplitude extracted from the data 
has been performed on the results marked "ISW-lensing sub- 
tracted". 10 Since the ISW-lensing bispectrum is peaked on 



10 More precisely, in the subtracted ISW-lensing results the equilateral 
and orthogonal primordial shapes are also fitted jointly, although this 
has a nearly negligible impact on the final result because the two shapes 
are by construction nearly perfectly uncorrected. 
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Table 9. Results for the /nl parameters of the primordial local, equilateral, and orthogonal shapes, determined by the KSW, binned 
and modal estimators from the SMICA, NILC, SEVEM, and C-R foreground-cleaned maps. Both independent single-shape results and 
results marginalized over the point source bispectrum and with the ISW-lensing bias subtracted are reported; error bars are 68% 
CL. 







Independent 




ISW-lensing subtracted 




JS.O W 


Binned 


iviouai 


KSW 


Binned 


Modal 


SMICA 














Local 


9.8 + 5.8 


9.2 ±5.9 


8.3 + 5.9 


2.7 + 5.8 


2.2 ± 5.9 


1.6 ±6.0 


Equilateral 


-37 + 75 


-20 ± 73 


-20 + 77 


-42 + 75 


-25 ± 73 


-20 ± 77 


Orthogonal 


-46 + 39 


-39 ± 41 


-36 + 41 


-25 + 39 


-17 ±41 


-14 ±42 


NILC 














Local 


11.6 + 5.8 


10.5 ± 5.8 


9.4 + 5.9 


4. S + 5 8 


j.kj m j .0 


2.7 + 6.0 


Equilateral 


-41 + 76 


-31 ±73 


-20 + 76 


-48 + 76 


-38 ± 73 


-20 ± 78 


Orthogonal 


HA -j- A(\ 
— / 4 + 4U 


— OZ ±41 


An -i- /in 
— OU ± 4U . . . 


-53 + 40 


-41 ±41 


-37 ± 43 


SEVEM 














Local 


10.5 + 5.9 


10.1 ±6.2 


9.4 + 6.0 


3.4 + 5.9 


3.2 ± 6.2 


2.6 ± 6.0 


Equilateral 


-32 + 76 


-21 ±73 


-13 + 77 


-36 + 76 


-25 ± 73 


-13 ±78 


Orthogonal 


-34 + 40 


-30 ± 42 


-24 + 42 


-14 + 40 


-9 ±42 


-2 ±42 


C-R 














Local 


12.4 + 6.0 


11.3 ±5.9 


10.9 + 5.9 


6.4 + 6.0 


5.5 ±5.9 


5.1 ±5.9 


Equilateral 


-60 + 79 


-52 ± 74 


-33 + 78 


-62 + 79 


-55 ± 74 


-32 ± 78 


Orthogonal 


-76 + 42 


-60 ± 42 


-63 + 42 


-57 + 42 


-41 ± 42 


-42 ± 42 



squeezed configurations, its impact is well known to be largest 
for the local shape. The ISW-lensing bias is also important for 
orthogonal measurements (there is a correlation coefficient r ~ 
-0.5 between the local and orthogonal CMB templates), while 
it is very small in the equilateral limit. The values of the ISW- 
lensing bias we subtract, summarized in Table 1, are calculated 
assuming the Planck best-fit cosmological model as our fidu- 
cial model. The same fiducial parameters were of course consis- 
tently used to compute the theoretical bispectrum templates and 
the estimator normalization. Regarding the point source contam- 
ination, we detect a Poissonian bispectrum at high significance 
in the SMICA map, see Sect. 5.3. However, marginalizing over 
point sources still carries a nearly negligible impact on the final 
primordial /nl results, because the Poisson bispectrum template 
has very small correlations with all the other shapes. 

In light of the discussion at the beginning of this section, we 
take the numbers from the KSW SMICA analysis in Table 8 as the 



Table 10. Results for the /nl parameters of the primordial local, 
equilateral, and orthogonal shapes, determined by the subopti- 
mal wavelet estimator from the SMICA foreground-cleaned map. 
Both independent single-shape results and results marginalized 
over the point source bispectrum and with the ISW-lensing bias 
subtracted are reported; error bars are 68% CL. As explained in 
the text, our current wavelets pipeline performs slightly worse in 
terms of error bars and correlation to primordial templates than 
the other bispectrum estimators, but it still provides a useful in- 
dependent cross-check of other techniques. 





Independent 


ISW-lensing subtracted 




Wavelets 


Wavelets 


SMICA 






Local 


10+8.5 


0.9 ± 8.5 


Equilateral . . . . 


89 + 84 


90 ± 84 


Orthogonal . . . . 


-73 + 52 


-45 ± 52 



final local, equilateral and orthogonal /nl constraints for the cur- 
rent Planck data release. These results clearly show that no evi- 
dence of NG of the local, equilateral or orthogonal type is found 
in the data. After ISW-lensing subtraction, all / N l for the three 
primordial shapes are consistent with at 68% CL. Note that 
these numbers have been cross-checked using two completely 
independent KSW pipelines, one of which is an extension to 
Planck resolution of the pipeline used for the WMAP analysis 
(Bennett et al. 2012). 

Unlike other methods, the KSW technique is not designed 
to provide a reconstruction of the full bispectrum of the data. 
However, the related skew-Q statistic described in Sect. 3.2.2 
allows, for each given shape, visualization and study of the con- 
tribution to the measured /nl from separate {-bins. This is a 
useful tool to study potential spurious NG contamination in the 
data. We show for the SMICA map in Fig. 5 the measured skew- 
Ce spectrum for optimal detection of primordial local, equilat- 
eral and orthogonal NG, along with the best-fitting estimates of 
/nl from the KSW method for different values of £, Contrary to 
the case of the point source and ISW-lensing foregrounds (see 
Sect. 5), the skew-Q statistics do not show convincing evidence 
for detection of the primordial shapes. In particular the skew- 
spectrum related to primordial local NG does not have the right 
shape, suggesting that whatever is causing this NG signal is not 
predominantly local. Again, point sources contribute very little 
to this statistic; ISW-lensing contributes, but only a small frac- 
tion of the amplitude, so there are indications of additional NG 
not captured by these foregrounds. In any event the estimators 
are consistent with no primordial signal of the types considered. 

As mentioned before, our analysis went beyond the simple 
application of the KSW estimator to the SMICA map. All /nl 
pipelines developed for Planck analysis were actually applied 
to all component-separated maps by SMICA, NILC, SEVEM, and 
C-R. We found from simulations in the previous Sections that 
the KSW, binned, and modal pipelines saturate the Cramer-Rao 
bound, while the wavelet estimator in its current implementation 
provides slightly suboptimal results. Wavelets remain however a 
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Fig. 5. Binned skew-Q statistics from the SMICA map for (a) 
local, (b) equilateral, and (c) orthogonal. Theoretical curves are 
not fitted to the data shown, but are plotted with the amplitude 
(the only free parameter) determined from the KSW technique. 
The dashed line shows the ISW-lensing contribution to the local 
statistic. There is no evidence for detection of primordial NG. 
Note that the error bars are underestimated, as they ignore data 
correlations. 



useful cross-check of the other methods, also given some techni- 
cal complementarities, e.g., they are the only approach that does 
not require inpainting, as explained in Sect. 3. Hence we include 
wavelet results, but only for SMICA. The / NL results for the opti- 
mal KSW, binned and modal bispectrum estimators, for the four 
component separation methods, are summarized in Table 9, one 
of the main products of our analysis of Planck data. The wavelet 
bispectrum analysis of SMICA is reported in Table 10. In the 
analysis, the KSW and binned bispectrum estimators considered 
multipoles up to { max = 2500, while the modal estimator went to 
Anax = 2000. As shown in Sect. 8. 1 and Table 16, error bars and 
central values for the three primordial shapes have converged at 
Anax = 2000, so the final primordial / NL estimates from the three 
pipelines are directly comparable. 11 

The binned bispectrum estimator used 51 bins, which were 
determined by optimizing the expected variance of the differ- 
ent /nl parameters, focusing in particular on the primordial 



11 The lower £ mslx for the modal pipeline is also a conservative choice 
in view of the large survey of "non-standard" models that will be pre- 
sented in Sect. 7.3 



shapes. The modal estimator employed a polynomial basis 
("max = 600) previously described in Fergusson et al. (2010a), 
but augmented with a local shape mode (approximating the 
SW large-angle local solution) to improve convergence in the 
squeezed limit. The above choices for the binned and modal 
methods produce a very high correlation (generally 99% or bet- 
ter) of the expanded/binned templates with the exact ones used 
by the KSW estimator. The wavelet estimator is based on third- 
order statistics generated by the different possible combinations 
of the wavelet coefficient maps of the SMHW evaluated at cer- 
tain angular scales. See for example Antoine & Vandergheynst 
(1998) and Martinez-Gonzalez et al. (2002) for detailed infor- 
mation about this wavelet. We considered a set of 15 scales 
logarithmically spaced between 1.3 and 956.3 arcmin and we 
also included the unconvolved map. The wavelet map w(Rj; b) 
(Eq. (60)) for each angular scale has an associated extended 
mask generated from the mask U73 following the procedure de- 
scribed and extensively used in Curto et al. 2009b, a, 2011a,b, 
2012; Donzelli et al. 2012; Regan et al. 2013. The wavelet coef- 
ficient maps are later combined into the third-order moments 
(Eq. (59)), for a total 816 different statistics, and these statistics 
are used to constrain / NL through a^ 2 test. 

The high level of agreement between results from the KSW, 
binned and modal /nl estimators, and from all the component 
separation pipelines, is representative of the robustness of our 
results with respect to residual foreground contamination, and 
is fully consistent with our preliminary MC analysis shown in 
Sect. 6. The scatter with wavelets is a bit larger, but this was 
expected due to the suboptimality of the wavelet estimator and 
is also in agreement with our MC expectations from the tests. 
Therefore wavelets do provide another successful cross-check. 

7.2. Bispectrum reconstruction 

As previously explained (see Sect. 3), in addition to looking in 
specific bispectrum-space directions and extracting the single 
number /nl for given shapes, the binned and modal pipelines 
have the capability to generate a smoothed (i.e., either coarse- 
grained in f-space, or truncated at a given expansion eigen- 
mode) reconstruction of the full bispectrum of the data. See also 
Planck Collaboration XXIII (2013). 

7.2.1. Modal bispectrum reconstruction 

The modal pipeline was applied to the Planck temperature 
maps for the foreground-separation techniques SMICA, NILC, 
and SEVEM (Fergusson et al. 2010a). For this analysis we used 
two alternative sets of hybrid basis functions in order to cross- 
check results and identify particular signals. First, we employed 
trigonometric functions (n m ax = 300) augmented with the SW 
local mode, together with the three separable modes contribut- 
ing to the CMB ISW-lensing signal. Secondly, we employed the 
same polynomial basis (n max = 600) with local SW mode as was 
used for / NL estimation. 

The modal coefficients extracted from the Planck SMICA, 
NILC, and SEVEM maps are shown in Fig. 8. Here we have used 
the hybrid Fourier modes with local and ISW-lensing modes. 



12 The boundary values of the bins are: 2, 4, 10, 18, 27, 39, 55, 75, 99, 
130, 170, 224, 264, 321, 335, 390, 420, 450, 518, 560, 615, 644, 670, 
700, 742, 800, 850, 909, 950, 979, 1005, 1050, 1110, 1150, 1200, 1230, 
1260, 1303, 1346, 1400, 1460, 1510, 1550, 1610, 1665, 1725, 1795, 
1871, 1955, 2091, 2240, and 2500 (i.e., the first bin is [2,3], the second 
[4,9], etc., while the last one is [2240,2500]). 
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Fig. 6. Full 3D CMB bispectrum recovered from the Planck foreground-cleaned maps, including SMICA (left), NILC (centre) and 
SEVEM (right), using the hybrid Fourier mode coefficients illustrated in Fig. 8, These are plotted in three-dimensions with multipole 
coordinates {t\, li, on the tetrahedral domain shown in Fig. 1 out to ^ max = 2000. Several density contours are plotted with red 
positive and blue negative. The bispectra extracted from the different foreground-separated maps appear to be almost indistinguish- 
able. 




Fig. 7. Planck CMB bispectrum detail in the signal-dominated regime showing a comparison between full 3D reconstruction using 
hybrid Fourier modes (left) and hybrid polynomials (right). Note the consistency of the main bispectrum properties which include 
an apparently 'oscillatory' central feature for low-^ together with a flattened signal beyond to t, < 1400. Note also the periodic CMB 
ISW-lensing signal in the squeezed limit along the edges of the tetrapyd. 



These amplitudes show remarkable consistency between the dif- 
ferent maps, demonstrating that the alternative foreground sepa- 
ration techniques do not appear to be introducing spurious NG. 
Note that here the coefficients are for the orthonormalized 
modes R„ (Eq. (63)) and they have a roughly constant variance, 
so anomalously large modes can be easily identified. It is ev- 
ident, for example, that among the low modes there are large 
signals, which include the ISW-lensing signal and point source 
contributions. 

Using the modal expansion of Eq. (45) with Eq. (63), we 
have reconstructed the full 3D Planck bispectrum. This is illus- 
trated in Fig. 6, where we show "tetrapyd" comparisons between 
different foreground cleaned maps. The tetrapyd (see Fig. 1) is 
the region defined by the multipoles that obey the triangle condi- 
tion, with £ < ( maK . The 3D plots show the reduced bispectrum of 
the map, divided by a Sachs-Wolfe CMB bispectrum solution for 



a constant primordial shape, S(k l ,k2,ki,) — I. This constant pri- 
mordial bispectrum template normalizaton is carried out in order 
to remove an ~ ( 4 scaling from the starting bispectrum (it is anal- 
ogous to multiplication of the power spectrum by ({( + 1)). To 
facilitate the interpretation of 3D bispectrum figures, note that 
squeezed configurations lie on the edges of the tetrapyd, flat- 
tened on the faces and equilateral in the interior, with bin on the 
diagonal. The colour levels are equally spaced with red denot- 
ing positive values, and blue denoting negative. Given the cor- 
respondence of the coefficients for SMICA, NILC, and SEVEM, 
the reconstructed 3D signals also appear remarkably consistent, 
showing similar contours out to I < 1500. At large multipoles ( 
approaching t max = 2000, there is increased randomness in the 
reconstruction due to the rise in experimental noise and some 
evidence for a residual point source contribution. 
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Fig. 8. Modal bispectrum coefficients for the mode expansion 
(Eq. (63)) obtained from Planck foreground-cleaned maps using 
hybrid Fourier modes. The different component separation meth- 
ods, SMICA, NILC and SEVEM exhibit remarkable agreement. The 
variance from 200 simulated noise maps was nearly constant for 
each of the 300 modes, with the average +l<x variation shown in 
red. 



o 



o 



0) 

a: 

o 





- NILC 




SEVEM 




ty - SMICA 














i JLlu i I i i_ 


_i i i i i i_ 





100 200 300 400 500 

Mode number n 



o 

CO 



o 
d 





- NILC 

SEVEM 
-SMICA 







200 400 600 800 1000 1200 1400 

Multipole I 

Fig. 9. The total integrated bispectrum F^ L defined in Eq. (64) 
as a cumulative sum over orthonormal modal coefficients /3* 2 
(upper panel) and over multipoles up to a given i (lower panel). 
Above, the relative quantity F^ h = F^ L - F^ h 2 is plotted, where 
F^ h 2 is the mean obtained from 200 CMB Gaussian maps with 
the standard deviation shown as the red line. Below the square 
of the bispectrum is integrated over the tetrapyd out to Z and its 
significance plotted relative to the Gaussian standard deviation 
( 1 cr red line) . A hybrid polynomial basis « max = 600 is employed 
in the signal-dominated region t < 1500. 



There are some striking features evident in the 3D bispec- 
trum reconstruction which appear in both Fourier and polyno- 
mial representations, as shown in more detail in Fig. 7. There is 
an apparent oscillation at low t < 500 already seen in WMAP-1 
(Fergus son et al. 2012). Beyond out to i ~ 1200 there are further 
distinct features (mostly "flattened" on the walls of the tetrapyd), 
and an oscillating ISW-lensing contribution can be discerned in 
the squeezed limit. Whatever its origin, Gaussian or otherwise, 
Fig. 7 reveals the CMB bispectrum of our Universe as observed 
by Planck. 

The cumulative sum F^ over the squared orthonormal co- 
efficients ffi 2 from Eq. (64) for the Planck data is illustrated in 
Fig. 9 (upper panel). The Planck bispectrum contribution can 
be directly compared with Gaussian expectations averaged from 
200 lensed Gaussian maps with simulated residual foregrounds. 
It is interesting to note that the integrated bispectrum signal 
fairly consistently exceeds the Gaussian mean by around 2cr 
over much of the domain. This includes the ISW and PS con- 
tributions for which subtraction only has a modest effect. Also 
shown (lower panel) is the corresponding cumulative F^ L quan- 
tity as a function of multipole i, for which features have visible 
counterparts at comparable ( in Fig. 7. Despite the high bispec- 
trum signal, this ^ 2 -test over the orthonormal mode coefficients 

is cumulatively consistent with Gaussianity. 

7.2.2. Binned bispectrum reconstruction 

As explained in Sect. 3.4.2, it is interesting to study the smoothed 
observed bispectrum divided by its expected standard devia- 
tion, since this will indicate if there is a significant deviation 
from Gaussianity for certain regions of ^-space. This quantity is 
shown in Figs. 10 and 1 1 as a function of (\ and (2, for two differ- 
ent values (or rather, bins) of (y. the intermediate value [610,654] 
in Fig. 10 and the high value [1330,1374] in Fig. 11. Each figure 
shows the results for the SMICA, NILC, SEVEM, and C-R cleaned 
maps as well as for the raw 143 GHz channel map. The bis- 
pectra were obtained with the binned bispectrum estimator and 
smoothed with a Gaussian kernel as explained in Sect. 3.4.2. 
Very blue or red regions indicate significant NG, regions that are 
less red or blue just represent expected fluctuations of a Gaussian 
distribution. 

From Fig. 1 at an intermediate value of {3 we can conclude 
that there is a very good agreement between SMICA, NILC, and 
SEVEM for all values of t\ and I2, and with C-R up to about 
l\,t2 ~ 1500. In fact, up to 1500 there is also a good agree- 
ment with the raw 143 GHz channel. We also see no significant 
non-Gaussian features in this figure (except maybe in the C-R 
and raw maps at (\ , (2 > 2000). 

Figure 1 1 at a high value of (3, on the contrary, shows signif- 
icant non-Gaussian features in the raw map, but much less NG in 
the cleaned maps. In particular one can see the point source bis- 
pectral signal at high-f approximately equilateral configurations. 
There is still an excellent agreement between SMICA, NILC, and 
SEVEM. The C-R map shows less NG than the other three cleaned 
maps, which is consistent with the absence of a detection of the 
Poisson point source bispectrum for C-R, see Table 3. 

7.3. Constraints on specific targeted shapes 

We have deployed the modal estimator to investigate a wide 
range of the inflationary models described in Sect. 2. This is 
the same validated estimator for which the standard /nl re- 
sults have been reported in the Sect. 7, but it is augmented with 
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Fig. 10. Smoothed observed bispectrum as determined with the binned estimator divided by its expected standard deviation, as a 
function of (\ and (2, with £3 in the bin [610,654]. From left to right on the top row are shown: SMICA, NILC, and SEVEM; and on the 
bottom row: C-R and the raw 143 GHz channel. 
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Fig. 11. Similar to Fig. 10, but with C 3 in the bin [1330,1374]. 



the primordial modal decomposition and projection described in 
Sect. 3.2.3. The resulting modal-projected local, equilateral and 
orthogonal shapes are ~99% correlated with those found using 



direct integration of Eq. (37) (as for the analysis above). Modal 
correlations for the other models investigated were determined 
for both the primordial shapes and the late-time projected de- 
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compositions and were all above 90%, unless stated otherwise. 
This primordial modal estimator pipeline has been applied al- 
ready extensively to the WMAP-7 data (Fergusson et al. 2012). 

7.3.1. Nonseparable single-field bispectrum shape results 

Having characterised single-field inflation bispectra using com- 
binations of the separable equilateral and orthogonal ansatze, 
we note that the actual leading-order non-separable contribu- 
tions (Eqs. (6, 7)) exhibit significant differences in the collinear 
(flattened) limit. For this reason we provide constraints on DBI 
inflation (Eq. (7)) and the two effective field theory shapes 
(Eqs. (5, 6)), as well as the ghost inflation bispectrum, which 
is an exemplar of higher-order derivative theories (specifically 
Eq. (3.8) in Arkani-Hamed et al. 2004). Using the primordial 
modal estimator, with the SMICA foreground-cleaned data, we 



find: 








/•DBI 

/nl _ 


: 11+69 


, F DBI-eq _ 
^ NL 


10 + 77) , 


f EFTl 


8+73 


z y^EFTl-eq 
( NL 


8+77), 


fEFTl 

/nl 


: 19 + 57 


/ ,-.EFT2-eq 
( NL 


27 ± 79) , 


^•Ghost 

/nl 


: -23 + 88 


/ nr^Ghost-eq 
^ NL 


-20 ± 75) . 



where we have normalized with the usual primordial /nl con- 
vention which is shape-dependent (i.e., the central value of the 
shape function is taken such that S(k,k,k) = 1). In parenthe- 
ses we also give a reweighted -F^l" 1 constraint for easier com- 
parison with the equilateral constraint from the same modal 
estimator, i.e., we have rescaled using the Fisher variance for 
the closely-related equilateral shape. Given the strong cross- 
correlation (above 95%) between all these models, the equi- 
lateral family results of (86) reveal larger differences around 
cr/3 than might be expected (and somewhat larger than ob- 
served previously in the WMAP data (Fergusson et al. 2012)). 
The reason for this variation between the equilateral shapes in 
Planck appears to be the additional signal observed in the flat- 
tened limit in the bispectrum reconstruction beyond the WMAP 
signal-dominated range (see Fig. 6). There is also a contribution 
from the small correlation difference between equilateral models 
from primordial modal and KSW methods. The results for these 
models for all the SMICA, NILC and SEVEM foreground-separated 
maps are given in Appendix B (Table B.3). 

7.3.2. Non-Bunch-Davies vacuum results 

We have investigated the non-separable shapes arising from ex- 
cited initial states (non-Bunch-Davies vacuum models) which 
usually peak in the flattened or collinear limit. In particular, we 
have searched for the four non-separable bispectra described in 
Eqs. (14) and (15), as well as the original flattened shape B^ BD 
(Eq. (6.2-3) in Chen et al. 2007b). This entails choosing suit- 
able cut-offs k c to ensure that the signal is strongly flattened 
(i.e., distinct from flat in Eq. (13)), while also accurately rep- 
resented by the modal expansion at both early and late times 
(Eqs. (54, 55)). For B^ BD , we adopted the same edge truncation 
and mild Gaussian filter described in Fergusson et al. (2012), 
while for B^bdi and b nbd2^ which are described by Eq _ (14 ) 5 

we chose k c = 0.001, and in Eq. (15) we take k c = 0.01. The 
shape correlations for most non-Bunch-Davies vacua were good 
(above 90%), except for the strongly squeezed model with os- 
cillations of Eq. (14) which was relatively poor (60%). Together 
with the orthogonal (Eq. (4)), flat (Eq. (13)) and vector (Eq. (19)) 
shapes, these non-Bunch-Davies models explore a broad range 



of flattened models, with a variety of different widths for picking 
out signals around the faces of the tetrapyd (see Fig. 1). 

The /nl results obtained for the non-Bunch-Davies models 
from the different foreground-cleaned map bispectra were con- 
sistent and the constraints from SMICA (for brevity) are given in 
Table 11. More comprehensive results from SMICA, NILC and 
SEVEM can be found in Table B.3 in Appendix B. Both B™ BD and 
^nbd2 produced raw results above 2<x, in part picking 

out the flattened signal observed in the bispectrum reconstruc- 
tion in Fig. 6. However, these flattened squeezed signals are also 
correlated with CMB ISW-lensing and so, after subtracting the 
predicted ISW bias (as well as the measured point source signal), 
most NBD /nl results were reduced to Icr or less (see "Clean 
/nl" column in Table 11). The exception was the most flattened 
model B^ BD which remained higher /^P B = 178 ± 78, i.e., with 
signals at 2.0cr, 1.8cr and 2.1cr for SMICA, NILC and SEVEM re- 
spectively. 

We emphasise that this has to be considered just as prelimi- 
nary study of flattened NG in the Planck data using four exem- 
plar models. In order to reach a complete statistical assessment 
of constraints regarding flattened models in forthcoming anal- 
yses, we will have to undertake a systematic search for best-fit 
Planck NBD models using the parameter freedom available. 

7.3.3. Scale-dependent feature and resonant model results 

We have investigated whether the Planck bispectrum reconstruc- 
tions include oscillations expected in feature or resonant models 
(Eqs. (16, 17)). Although poorly correlated with scale-invariant 
shapes, the feature and resonant models have (at least) two free 
parameters - the period k c and the phase <p - forming a model 
space which must be scanned to determine if there is any sig- 
nificant correlation (in the absence of any physical motivation 
for restricting attention to specific periodicities). We have under- 
taken an initial survey of these models with the wavelength range 
defined by the native resolution of the present modal estimator 
(hybrid local polynomials with 600 modes), similar to the feature 
model search in WMAP data in Fergusson et al. (2012). For fea- 
ture models of Eq. (16) we can obtain high correlations (above 
95%) for the predicted CMB bispectrum if we take k c > 0.01, 
that is, for an effective multipole periodicity ( c > 140 feature 
models are accurately represented. 

The results of a first survey of feature models in the Planck 
data is shown in Table 12 for 0.01 < k c < 0.1 and phases 
4> = 0, 7r/4,7r/2,37r/4 (for <\> > n we will identify a correlation 
with the opposite sign). Again, there was good consistency be- 
tween the different foreground-separation methods SMICA, NILC 
and SEVEM showing that the results are robust to potential resid- 
ual foreground contamination in the data. For brevity we only 
give SMICA results here, while providing measurements from 
other component separation methods in Appendix B. Feature 
signals are typically largely uncorrelated with the ISW-lensing 
or point sources, but nevertheless we subtract these signals and 
give results for the cleaned /nl- The Table 12 results show that 
there is a parameter region around 0.01 < k c < 0.025 for which 
signals well in excess of 2<x are possible (we undertook a broader 
search with 0.01 < A: c < 0.1 but found only a low signal be- 
yond k > 0.3). It appears that some feature models are able to 
match the low-^ 'plus-minus' and other features in the Planck 
bispectrum reconstruction (see Fig. 6). The best fit model has 
k c = 0.0185 (£ c « 260) and phase = with a signal -3cr. 
As a further validation step of our results, we also re-analysed 
the models with > 2.5cr significance using a different modal de- 
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Table 11. Constraints on flattened or collinear bispectrum models (and related models) using the SMICA foreground-cleaned Planck 
map. These bispectrum shapes, with equation numbers given, are described in detail in the text. 



Flattened model (Eq. number) 


Raw /nl 


Clean / NL 


A/nl 


a 


Clean <x 


Flat model (13) 


70 


37 


77 


0.9 


0.5 


Non-Bunch-Davies (NBD) 


178 


155 


78 


2.2 


2.0 


Single-field NBD1 flattened (14) 


31 


19 


13 


2.4 


1.4 


Single-field NBD2 squeezed (14) 


0.8 


0.2 


0.4 


1.8 


0.5 


Non-canonical NBD3 (15) 


13 


9.6 


9.7 


1.3 


1.0 


Vector model L = 1 (19) 


-18 


-4.6 


47 


-0.4 


-0.1 


Vector model L = 2 (19) 


2.8 


-0.4 


2.9 


1.0 


-0.1 



Table 12. Planck bispectrum estimation results for feature models compared to the SMICA foreground-cleaned maps. This prelim- 
inary survey on a coarse grid in the range 0.01 < k c < 0.025 and < k c < n/4 finds specific models with significance up to 
99.7%. 



Phase = = tt/4 <p = n/2 <f> = 3tt/4 



Wavenumber 


,/.NL 


± A/nl 


(o-) 


/nl 


± A/nl 


(o-) 


./.NL 


± A/nl 


(o-) 


/nl 


± A/nl 


(o-) 


fe c 


= 0.01000 


-110 + 


159 (- 


■0.7) 


-98 ± 


167 (- 


-0.6) 


-17 ± 


147 (- 


-0.1) 


56 ± 


142 ( 


0.4) 


fee 


= 0.01125 


434 + 


170 ( 


2.6) 


363 ± 


185 ( 


2.0) 


57 ± 


183 ( 


0.3) 


-262 ± 


168 (- 


1.6) 


fee 


= 0.01250 


-70 + 


158 (- 


•0.4) 


130 ± 166 ( 


0.8) 


261 ± 


167 ( 


1.6) 


233 + 


159 ( 


1.5) 


fee 


= 0.01375 


35 + 


162 ( 


0.2) 


291 ± 


145 ( 


2.0) 


345 ± 


147 ( 


2.3) 


235 ± 


162 ( 


1.5) 


fee 


= 0.01500 


-313 + 


144 (- 


-2.2) 


-270 ± 


137 (- 


-2.0) 


-95 ± 


145 (- 


-0.7) 


179 + 


154 ( 


1.2) 


fee 


= 0.01625 


81 + 


126 ( 


0.6) 


177 + 


141 ( 


1.2) 


165 ± 


144 ( 


1.1) 


51 ± 


129 ( 


0.4) 


fee 


= 0.01750 


-335 + 


137 (- 


-2.4) 


-104 + 


128 (- 


0.8) 


181 ± 


117 ( 


1.5) 


366 ± 


126 ( 


2.9) 


fee 


= 0.01875 


-348 + 


118 (- 


-3.0) 


-323 ± 


120 (- 


2.7) 


-126 ± 


119 (- 


-1.1) 


137 ± 


117 ( 


1.2) 


fee 


= 0.02000 


-155 + 


110 (- 


-1.4) 


-298 ± 


119 (- 


-2.5) 


-241 ± 


113 (- 


■2.1) 


-44 ± 


105 (- 


0.4) 


fee 


= 0.02125 


-43 


+ 96 (- 


-0.4) 


-186 ± 


107 (- 


1.7) 


-229 + 


115 (- 


-2.0) 


-125 + 


104 (- 


1.2) 


fee 


= 0.02250 


22: 


t95 ( 


0.2) 


-115 


+ 92 (- 


1.2) 


-194 ± 


105 (- 


-1.8) 


-148 ± 


107 (- 


1.4) 


fee 


= 0.02375 


70 + 


100 ( 


0.7) 


-56 


+ 94 (- 


0.6) 


-159 


±93 (- 


-1.7) 


-164 ± 


101 (- 


1.6) 


fee 


= 0.02500 


106 : 


t93 ( 


1.1) 


6: 


t97 ( 


0.1) 


-103 


+ 98 (- 


■1.1) 


-153 


+ 94 (- 


1.6) 



composition, namely an oscillating Fourier basis (n max = 300) 
augmented with a local SW mode (the same used for the recon- 
struction plots in Sect. 7). The results from this basis are shown 
in Appendix B and they are fully consistent with the polynomial 
measurements presented here. The previous best-fit WMAP fea- 
ture model, k c = 0.014 (i c ss 200) and phase <f> = 3n/4, attained 
a 2.15cr signal with £ < 500 (Fergusson et al. 2012), but it only 
remains at this level for Planck. 

We note however that the apparently high statistical signifi- 
cance of these results is much lower if we consider this to be a 
blind survey of feature models, because we are seeking several 
uncorrelated models simultaneously. Following what we did for 
our study of impact of foregrounds in Sect. 8, we considered a 
set of 200 realistic lensed FFP6 simulations, processed through 
the SMICA pipeline, and including realistic foreground residuals. 
If we use this accurate MC sample to search for the same grid 
of 52 feature models as in Table 12, we find a typical maximum 
signal of 2.23(+0.56) cr. Searching across all feature models (see 
below) studied here yields an expected maximum 2.37(±0.53) <x 
(whereas the survey for all 511 models from all paradigms in- 
vestigated yielded 2.55(±0.52) cr). This means that our best-fit 
model from data has a statistical significance below 1 .5 cr above 
the maximum signal expectation from simulations, so we con- 
clude that we have no significant detection of feature models 
from Planck data. 

Feature models typically have a damping envelope repre- 
senting the decay of the oscillations as the inflaton returns to 
its background slow-roll evolution. Indeed, the feature envelope 
is a characteristic of the primordial mechanism producing the 
fluctuations, decaying as k increases for inflation while rising 




Fig. 12. CMB bispectrum shown for the best-fit feature model 
with an envelope with parameters k = 0.01875, phase = and 
Ak = 0.045 (see Table 13). Compare with the Planck bispectrum 
reconstruction, Fig. 7. 

for contracting models like the ekpyrotic case (Chen 2011). We 
have made an initial survey to determine whether a decaying 
envelope improves the significance of any feature models. The 
envelope employed was a Gaussian centred at k c = 0.045 with 
a falloff Ak = 0.015, 0.02, 0.025, 0.03, 0.035, 0.04, 0.045 and re- 
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Table 13. Feature model results with an envelope decay function. Results are only presented for feature models with better than 
95% CL result on the full domain (see Table 12). 
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fee 
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-665 ± 
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2.9) 


-593 ± 


185 (- 


3.2) 


-500 + 


160 (- 


3.1) 


-298 + 


119 (- 


•2.5) 



suits for specific parameters are given in Table 13. The best fit 
model remains k = 0.01875 {£ c = 265) with phase cp — and 
the significance rises to 3.23<x, together with a second model 
k = 0.02 (t c = 285) cp = n y '4. However the caveats about blind 
survey statistics previously noted also do not allow a claim of 
any detection in this case. A plot of the best-fit feature model 
with a decay envelope is shown in Fig. 12, for which the main 
features should be compared with those in Fig. 7. Non-Gaussian 
bispectrum signals from feature models typically produce coun- 
terparts in the power spectrum as will be described in Sect. 9. An 
improved statistical interpretation of the results presented in this 
Section will be possible when this additional investigation will 
be completed. 

We have also undertaken a survey of resonant models and 
the non-Bunch-Davies resonant models (or enfolded resonance 
models). With the modal estimator, we can achieve high ac- 
curacy for the predicted bispectrum for k c > 0.001 (note that 
this has a different logarithmic dependence to feature models 
and a varying effective i c ). For the resonance model shape of 
Eq. (18), we have not undertaken an extensive survey, except 
selecting a likely range for a high signal with periodicity com- 
parable to the feature model, that is, with 0.25 < k c < 0.5 
and phases = 0, n/4, k/2, 3n/4, n. However, no signif- 
icant signal was found (all below lcr), as can be verified in 
Table B.l in Appendix B. For the enfolded resonance model 
shape of Eq. (18) , we have undertaken a preliminary search in 
the range 4 < k c < 12 with the same phases. Again, no signifi- 
cant signal emerges from the Planck data, as shown in Table B.2 
in Appendix B. 

7.3.4. Directional dependence motivated by gauge fields 

We have investigated whether there is significant NG from 
bispectrum shapes with non-trivial directional dependence 
(Eq. (19)), which are motivated by inflationary models with vec- 
tor fields. Using the primordial modal estimator we obtained a 
good correlation with the L = 1 flattened type model, but the 
squeezed L = 2 model produced a relatively poor correlation 
of only 60%, given the complexity of the dominant squeezed 
limit. Preliminary constraints on these models are given in the 
Table 11, showing no evidence of a significant signal. 

7.3.5. Warm inflation 

Warm inflation produces a related shape with a sign change 
in the squeezed limit. This also had a poor correlation, until 
smoothing (WarmS) was applied as described in Fergusson et al. 
(2012). The resulting bispectrum shows no evidence for signifi- 
cant correlation with Planck data (SMICA), 

^armS = 4 + 33 . (87) 



The full list of constraints for SMICA, NILC and SEVEM models 
can be found for warm inflation and vector models in Table B.3 
in Appendix B . 

7.3.6. Quasi-single-field inflation 

Finally, quasi-single-field inflation has been analysed constrain- 
ing the bispectrum shape of Q (Eq. (12)), that depends on two 

osi 

parameters, v and f^ L . In order to constrain this model we have 

calculated modal coefficients for < v < 1.5 in steps of 0.01 

(so 151 models in total). These were then applied to the data 

and the one with the greatest significance was selected. Results 

are shown in Fig. 24. The maximum signal occurred at v = 1.5, 
osi 

f^i = 4.79 (0.3 lcr). To obtain error curves we performed a 
full likelihood using 2 billion simulations following the method 
described in Sefusatti et al. (2012). Such a large number of sim- 
ulations was possible as they were generated from the modal (3- 
covariance matrix which is calculated once from the 200 Planck 
realistic CMB simulations, rather than repeatedly from the CMB 
simulations themselves. The procedure is to take the 151 X 151 
correlation matrix for the models (this is just the normalized dot 
product of the modal coefficients). This is then diagonalised us- 
ing PCA, after which only the first 5 eigenvalues are kept as the 
remaining eigenvalues are < 10~ 10 . The y6-covariance matrix is 
projected into the same sub-basis where it is also diagonalised 
via PCA into 5 orthonormal modes, with the two leading modes 
closely correlated with local and equilateral. The procedure by 
which to produce a simulation is to generate five Gaussian ran- 
dom numbers and add the mean values obtained from the Planck 
data, rotating them to the sub-basis where we determine the v 
with the greatest significance. The result is then projected back 
to the original space to determine the related /nl- The two billion 
results from this MC analysis are then converted into confidence 
curves plotted in Fig. 24. The curve shows that there is no pre- 
ferred value for v with all values allowed at 3cr. This reflects the 
results obtained from data previously, where we found the least 
preferred value of v = 0.86 had only a marginally lower signif- 
icance of 0.28<x (Sefusatti et al. 2012). Of course, these conclu- 
sions are directly related to the null results for both local and 
equilateral templates. 

7.4. Constraints on local non-Gaussianity with Minkowski 
Functionals 

In this Subsection, we present constraints on local NG ob- 
tained with Minkowski Functionals (MFs). MFs describe the 
morphological properties of the CMB field and can be used as 
generic estimators of NG (Komatsu et al. 2003; Eriksen et al. 
2004; De Troia et al. 2007; Hikage et al. 2008; Curto et al. 2008; 
Natolietal. 2010; Hikage & Matsubara 2012; Modest et al. 
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Fig. 13. The Wiener filter W M used to constrain f^ di with MFs. 

Table 14. Validation tests with MFs: results for /^ al obtained 
using the filter W M , for £ max = 2000 and V side = 2048. 
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0.21 d 


10.71 


Realistic Noise, = 12 


12.14 d 


- 13.12 


Mask (/ sty = 0.73), realistic noise / t | 1 ° L cal = 12 


12.18 d 


.- 18.13 



7.4.1. Method 

We constrain /}Jj£ using the optimized procedure described in 
Ducout et al. (2013). To obtain constraints on f^f 1 , we apply a 
specific Wiener filter on the map (Wm), shown in Fig. 13. We do 
not use here the filter designed to enhance the information from 
the gradients of the map {Wr>i = + 1)Wm), because this 

filter is very sensitive to small-scale effects and may be biased 
by foreground residuals. 

We use maps at HEALPix resolution Af s jd e = 
1024 (Gorskietal. 2005) and f m . dx = 2000. Our results 
are based on the four normalized 13 functionals v k (k = 0,3) 
(respectively Area, Perimeter, Genus and V c i ust er), computed on 
n t h = 26 thresholds v, between v m ; n = -3.5 and v max = +3.5 in 
units of the standard deviation of the map. 

We combine all functionals into one vector y (of size n - 
104). We then analyse this vector in a Bayesian way to obtain 
a posterior for the , and hence an estimate of this parame- 
ter. The principle is to compare the vector measured on the data 
y to the ones measured on non-Gaussian simulations with the 
same systematic effects (realistic noise, effective beam) and data 
processing (Wiener filtering) as the data, yif^ 31 )- Modelling the 
MFs as multivariate Gaussians we obtain the posterior distribu- 
tion for /'° cal with ax 2 test : 



P(f^\y)< 



exp 



(88) 



2013). As they are sensitive to every order of NG, they 
can be used to constrain different bispectrum and trispectrum 
shapes (Hikage et al. 2006, 2008; Hikage & Matsubara 2012). 
They are therefore complementary to, and a useful valida- 
tion of, optimal estimators. Their precise definition and ana- 
lytic formulations are presented in Planck Collaboration XXIII 
(2013). The MF technique is also used in the companion paper 
Planck Collaboration XXV (2013). 

We review here the properties of MFs, as a complementary 
tool to poly-spectrum based estimators. 

First, they are defined in real space, which makes MFs ro- 
bust to masking effects and no linear term is needed to take 
into account the anisotropy of the data model. Second, as MFs 
are sensitive to every non-Gaussian feature in the maps, they 
can be a useful probe of every potential bias in the bispectrum 
measurement, in particular the different astrophysical contami- 
nations (foregrounds and secondaries). 

There is a limitation to MF studies: they can be expressed in 
terms of weighted sums of the bispectrum (and trispectrum) in 
harmonic space (Matsubara 2010), hence the angle-dependence 
of the bispectrum is partially lost. This makes MFs subopti- 
mal in two ways: increasing error bars for constraints on spe- 
cific shapes and reducing the distinguishability of different bis- 
pectrum shapes. This lack of specificity can introduce biases, 
as MFs will partially confuse primordial and non-primordial 
sources of NG and can introduce degeneracies between differ- 
ent primordial shapes. Constraints on orthogonal and equilateral 
shapes are quite degenerate with MFs, we therefore chose here to 
focus on the local bispectrum shape. We also leave trispectrum 
analyses for future studies. 

An attractive feature of MFs is their linearity for weak 
NG (/nl) and weak signals (such as point sources, and 
Galactic residuals after masking and component separation) 
(Ducout et al. 2013). This property can be used to estimate dif- 
ferent known non-primordial contributions. 



with 



y 



y 



■yifv? 1 ) 



(89) 



Since NG is weak, the covariance matrix C is computed with 
10 4 Gaussian simulations, again reproducing effective beam, re- 
alistic noise and filtering of the data. The dependence of the MFs 
on /nl ' ^(/nT 3 '-*' * s obtained as an average of y measured on 100 
simulations. The simulations used here are based on the WMAP- 
7 best-fit power spectrum (Komatsu et al. 2011), using the pro- 
cedure described in Eisner & Wandelt (2009). 



7.4.2. Validation tests 



We report here validation of the MFs estimator on in 
thoroughly realistic Planck simulations. This validation subsec- 
tion is analogous to Sect. 6. 1 concerning bispectrum-based es- 
timators. The same tests (ideal Gaussian maps, full-sky non- 
Gaussian maps with noise and non-Gaussian maps with noise 
and mask) are performed, but different non-Gaussian simula- 
tions are used. Non-Gaussian CMB simulations as processed in 
Fergusson et al. (2010a) only guarantee the correctness of the 3- 
point correlations. Since the MFs are sensitive to higher-order 
n-point functions, they were validated with physical simulations 
(Eisner & Wandelt 2009). 

The first test consists of 100 simulations of a full-sky 
Gaussian CMB, with a Gaussian beam with FWHM of 5 arcmin 
and without any noise, cut off at { max = 2000, with N s id e = 2048. 
Here validation tests were made at N s a e = 2048, but results (es- 
timate and error bars) remain the same at N s ife = 1024 as we 
keep the same ^ max . The second test includes non-Gaussian sim- 
ulations with = 12 and realistic coloured and anisotropic 



13 Raw MFs VV depend on the Gaussian part of fields through a nor- 
malization factor At that is a function only of the power spectrum shape. 
We therefore normalize functionals v k = VklA k to focus on NG, see 
Planck Collaboration XXIII (2013) and references therein. 



33 



Planck Collaboration: Planck 2013 Results. XXIV. Constraints on primordial NG 



Oo 

I 

o 




SMICA 




0.73, W M 



On 
I 




4 ♦ SMICA 


i 1 i i i 1 ■■ i 


™ flocal FTfl 

JNL — 0U 


; 








\/ / sky = 0.73, W u J 

l l 


" (d) 


-2 


2 



Fig. 14. The MFs curves for SMICA at jV side = 1024 and { max = 2000, for the four functionals u k : (a) Area, (b) Perimeter, (c) 
Genus, and (d) N c \ uster . The curves are the difference of each normalized MF, measured from the data, to the average from Gaussian 
Planck realistic simulations (not lensed). The difference curves are normalized by the maximum of the Gaussian curve. To compare 
the curves to the presence of primordial NG, the average difference curves for non-Gaussian simulations with f^ dl = 50 is also 
represented (100 simulations). 



Table 15. Estimates of obtained with MFs on Planck data. Foreground and secondary effects are evaluated in terms of /J 
Results are for SMICA at N side = 1024 and ? m . dx = 2000. 



local 





•'NL 


ai 


Source 


Corresponding A/jJj]" 1 


Raw map 


19.1 + 


19.3 






Lensing subtracted 


8.5 ± 


20.5 


Lensing 


+ 10.6 


Lensing+PS subtracted 


7.7 ± 


20.3 


Point sources 


+0.8 


Lensing+CIB subtracted 


7.5 ± 


20.5 


CIB 


+ 1.0 


Lensing+SZ subtracted 


6.0 ± 


20.4 


sz 


+2.5 


All subtracted 


4.2 ± 


20.5 


All 


+ 14.9 



noise, processed through the Planck simulation pipeline and the 
component-separation method SMICA. Finally, in the third test 
we add the union mask U73 to the previous simulations, mask- 
ing both the Galaxy and the brightest point sources, and leaving 
73% of the sky unmasked. Only the inpainting of the smallest 
holes in the point sources part of the mask was performed. For 
these three tests, the results are presented in Table 14. We give 
here the average estimate and error bar obtained on the 100 simu- 
lations, when we use the four functionals. The results show that 
the MF estimator is unbiased, robust, and a competitive alter- 
native to bispectmm-based estimators. Moreover a map-by-map 



comparison of the results obtained on ffi£ with KSW and MFs 
estimators showed a fair agreement between the two methods. 

7.4.3. Results 

For our analysis we considered a foreground-cleaned map 
obtained with the component separation method SMICA at 
Nside = 1024 and ^ max = 2000. As for the previous results 
in this Section, we used the union mask U73, which leaves 
/sky = 73% of the sky after masking Galaxy and point sources. 
To take into account some instrumental effects (asymmetry 
of beams, component separation processing) and known non- 
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Gaussian contributions (lensing), we used realistic unlensed and 
lensed simulations (10 3 ) of Planck data (FFP6 simulations, see 
Planck Collaboration ES 2013). First, MFs were applied to the 
unlensed simulations and the resulting curves served to calibrate 
the estimator, as the Gaussian part of the NG curves yif^ 1 ) 14 . 
This estimate is referred to as the "raw map". Secondly, MFs 
were applied to the lensed simulations, and the same procedure 
was applied, the result being referred as "lensing-subtracted". 
We summarize the procedure in the following equation: 



</(/NL al 



^ = w G + Aw NG 

' "Planck simulations, lensed "/nli NG simulations 



(90) 



Here we assume that the MF respond linearly to lensing at first 
order and that primordial NG and lensing contributions are there- 
fore additive. 

Additionally, we tried to characterize other non-primordial 
contributions that one could expect in masked SMICA-cleaned 
maps covering 73% of the sky. To this end, we used simula- 
tions of extragalactic foregrounds and secondary anisotropics: 
uncorrelated (Poissonian) point sources; clustered CIB; and 
SZ clusters. These component simulations reproduce accurately 
the whole Planck data processing pipeline (beam asymme- 
try, component separation method). Using the linearity of MFs 
(Ducout et al. 2013), we could introduce these effects as a simple 
additive bias on the curves following 



y 



^subtracted + A -PS + A ~C1B + A ~SZ 



(91) 



where Aj/ PS '" is the average bias measured on 100 simulations. 
Note that the SZ simulation does not take into account the SZ- 
lensing correlation, which is expected to be negligible given the 
error bars. 

Results are summarized in Table 15 and MFs curves are 
shown in Fig. 14, without including the lensing subtraction 
("raw curves"). Considering the larger error bars of MF estima- 
tors, the constraints obtained are consistent with those from the 
bispectrum-based estimators, even without subtracting the ex- 
pected non-primordial contributions. Moreover, results are quite 
robust to Galactic residuals: constraints obtained with other 
component separation methods (NILC and SEVEM), with differ- 
ent sky coverage, differ from the SMICA results presented here 
by less than A/^ al = 1. 

7.5. CMB trispectrum results 

As shown in Fig. 15 the modulation reconstruction mean field 
has two large contributions, one from the mask and one from 
anisotropic noise, reflecting the fact that they both look like a 
large spatially-varying modulation of the fluctuation power. The 
noise we use to estimate the mean field is taken from FFP6 sim- 
ulations, adjusted with an additional lOyuKarcmin white noise 
component to match the power spectrum in the observed maps. 
However this is still only an approximate description of the in- 
strumental noise present in the data. The mean field from non- 
independent maps (e.g., 143 x 143 GHz maps) shows a large 
noise anisotropy that is primarily quadrupolar before masking, 
and any mismatch between the simulated noise and reality would 
lead to a large error in the mean field subtraction. By instead us- 
ing the modulation estimator for 143 x 2 17 GHz maps errors due 
to misestimation of the noise are avoided, and the mean field 
is then dominated by the shape of the Galactic cut, which is 
well known, and a smaller uncertainty from assumed simulation 




3 4 



10 
L 



100 



14 The overall effect of data processing on the /jjjj al constraint from 
MFs was evaluated as /^'"'(process.) ~ 3. 



Fig. 15. The two upper maps show the modulation reconstruc- 
tion mean field / MF (n) at L < 100 , which is essentially a map 
of the expected total small-scale power on the masked map as 
a function of position (assuming there is no primordial power 
modulation). The top mean field map from the 143 GHz auto 
estimator has a large signal from both the cut (which can be 
calculated accurately), and from the noise anisotropy (aligned 
roughly with the ecliptic, which cannot be estimated very accu- 
rately from simulations). The lower mean field is the 143 x 217 
GHz cross-estimator map, and does not have the contribution 
from the noise anisotropy (note the colour scale is adjusted). 
The lower plot shows the corresponding mean field power spec- 
tra compared to the reconstruction noise (connected part of 
the trispectrum); the reconstruction noise is much smaller than 
both the detector noise and mask contributions to the mean field. 
Since any tnl signal is all on large scales we do not reconstruct 
modes above L max = 100. 



power spectrum and beam errors (see Fig. 15). For this reason for 
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Fig. 16. Power spectrum of the power modulation reconstructed 
from 143 x 217 GHz maps. Shading shows the 68%, 95% and 
99% CL intervals from simulations with no modulation or kine- 
matic signal. The dashed lines are when the mean field simula- 
tions include no kinematic effects, showing a clear detection of a 
modulation dipole. The blue points show the expected kinematic 
modulation dipole signal from simulations, along with l<x error 
bars (only first four points shown for clarity). The solid line sub- 
tracts the dipolar kinematic signal in the mean fields from sim- 
ulations including the expected signal, and represents out best 
estimate of the non-kinematic signal (note this is not just a sub- 
traction of the power spectra since the mean field takes out the 
fixed dipole anisotropy in real space before calculating the re- 
maining modulation power). The dotted line shows the expected 
signal for tnl = 1000. 



our main result we work with modulation reconstructions gener- 
ated from 143 x 217 GHz maps with independent noise, which 
removes the leading error due to noise mean field misestimation. 

Figure 16 shows the reconstructed modulation power from 
143 x 217 GHz maps that we use for our analysis. We show 
two results: one where we do not include the expected kinematic 
dipole signal in the mean field subtraction (see Sect. 5.4), and 
one were we do so that the reconstruction should then be domi- 
nated by the primordial and any unmodelled systematic effects. 
In the first case the 143 x 217 result gives a clear first detection 
of the dipolar kinematic modulation signal of roughly the ex- 
pected magnitude (see Planck Collaboration XXVII (2013) for a 
more detailed discussion of this signal). Including the expected 
kinematic signal in the simulations (and hence the mean field) 
removes this signal, giving a cosmological modulation recon- 
struction that is broadly consistent with no modulation (statisti- 
cal isotropy) except for the anomalous very significant signal in 
the modulation octopole. 

Note that only the two-point reconstruction is free from noise 
bias, the four-point is still sensitive to noise modelling at the 
level of the subtraction of the A^ 0) (Eq. (72)) reconstruction noise 

power spectrum. However as shown in the Fig. 16, A^ 0) is not 
that much larger than the reconstruction scatter at low multipoles 
where the tnl signal peaks, so the sensitivity to noise misestima- 
tion is much less than in the mean field subtraction (where the 




2000 



Fig. 17. Distribution of f NL estimators from Gaussian simula- 
tions (L max = 10) compared to data estimates (vertical lines). 
The distribution is rather skewed because the main contributions 
are from L < 4 where the modulation power spectra have skewed 
X 2 distributions with low degrees of freedom. The red line shows 
the predicted distribution for a weighted sum of tnl(^) estima- 
tors assuming the reconstructed modulation modes are Gaussian 
with 2L + 1 modes measured per L, which fits the full simula- 
tions well. The black vertical lines show the data estimates from 
Lmax = 10, and should be compared to the green which have 
Lmax = 2 and hence are insensitive to the anomalous octopole 
signal. The dashed lines are Tnl,i, the slightly more optimal vari- 
ant of the estimator. 



large-scale noise anisotropy gives a large-scale mean field in the 
auto-estimators orders of magnitude larger, Fig. 15). 

The t nl estimator from the 143 X 217 GHz modulation re- 
construction gives tnl = 442, compared to a null hypothesis 
distribution -452 < t nl < 835 at 95% CL (cr TNL * 335). Our 
quoted error bar is assuming zero signal so that there is no signal 
cosmic variance contribution, and the bulk of the apparent signal 
is coming from the high octopole seen in Fig. 16. The alterna- 
tive estimator Tnl,i gives a slightly different weighting to the 
octopole, giving tnl.i = 569 with an expected null-hypothesis 
<x TNL ~ 332. The surprisingly large difference between the es- 
timators can be explained as due to the large octopole signal, 
which has tnl(£ = 3) ~ 6000. However the shape of the total 
signal would not be expected from a genuine tnl signal, since 
as shown in Fig. 16 on average the latter is expected to fall off 
approximately proportional to l/L 2 (i.e., a large primordial tnl 
would be expected in most realisations to give large dipole and 
quadrupole signals that we do not see). If we estimate r NL using 
Lmax = 2 we obtain tnl = 165 with only a slightly larger null- 
hypothesis error cr TNL = 349, where in this case tnl,i = 172. 

Note that the distribution of tnl is quite skewed because 
the signal is dominated by the very-low multipole modulation 
power spectra which have skewed ^ 2 -like distributions due to 
the large cosmic variance there (see Fig. 17; Hanson & Lewis 
2009; Smith & Kamionkowski 2012). The reconstruction noise 
acts like nearly-uncorrelated Gaussian white noise, so each of 
the c£ comes from a sum of squares of ~ 2L + 1 modulation re- 
construction modes; the shape of the ?nl distribution is consis- 
tent with what would come from calculating a weighted sum of 
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Fig. 18. Approximate posterior distributions P(tnl\tnl(L)) for a 
range of L max . The distributions have broad tails to high values 
because of the small number of large-scale modulation modes 
that are measured, and hence large cosmic variance. For L max > 
3 the distributions are pulled away from zero by the significant 
octopole modulation signal observed, and are gradually move 
back towards zero as we include more modulation modes that 
are inconsistent with large tnl values. As shown in Fig. 19 the 
octopole has significant frequency dependence and is therefore 
unlikely to be physical. 

^-distributed random variables. If we assume that any primor- 
dial modulation modes giving rise to a physical tnl signal are 
also Gaussian, for any given physical tnl the tnl(^) estimators 
would also have x 2 distributions. This allows us to evaluate the 
posterior distribution of tnl given the observed tnl, in exactly 
the same way that one can do for the CMB temperature power 
spectrum. For each L the posterior distribution P(tnl(L)\tnl(L)) 
on the full sky would have an inverse-gamma distribution. We 
follow Hamimeche & Lewis (2008) by generalizing this to a cut- 
sky approximation for a range of multipoles: 

-21nP({r NL (L)}|{r NL (L)})* 

J] [g(x(L))N<®(L)] [M-\ L , [ti™ (L')g(x(L'))] (92) 

LL' 

where Miu is the covariance of the estimators calculated from 
null hypothesis simulations, Nt°J l (L) = kjjslf /C^ is the estima- 
tor reconstruction noise, 



x(L) = 



f NL (L) +<[(£) 



(93) 



and 



g(x) = sign(* - 1) yj2(x - ln(x) - 1). 



(94) 

For uncorrected ^-distributed estimators this distribution re- 
duces to the exact distribution, a product of inverse-gamma dis- 
tributions. This approximation to the posterior can be used to 
evaluate the probability of any scale dependent r NL (L), and does 
not rely on compression into a single tnl estimator; it can there- 
fore fully account for the observed distribution of modulation 
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Fig. 19. Comparison of the un-normalized modulation power 
C L with various combinations of frequencies. The middle plot 
shows the results used for our main Tnl result, since it removes 
the significant largely-quadrupolar signal from anisotropic noise 
misestimation seen in the two other plots. The noticeable differ- 
ence in the odd octopole signal between channels indicates that 
the residual signal in 143 x 217 GHz is unlikely to be physical, 
but we cannot currently identify its origin. 



power between L. Here we focus on the main case of interest 
where tnl(^) is nearly-scale-invariant so that for all L we have 
?nl(£) = tnl- 

The resulting posterior distributions of tnl are shown in 
Fig. 18 for a range of L max . These are strongly skewed, in the 
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same way as the posterior from the low quadrupole in the CMB 
temperature data. The high octopole is pulling the distributions 
up to higher tnl values, but increasing L max can reduce the high- 
ly tail because very large tnl values are inconsistent with the low 
modulation power seen at L + 3. With L max = 2 the posterior 
peaks near zero, but the distribution is then very broad because 
there are only about 8 modes, which therefore have large cosmic 
variance. For L max = 50 we find tnl < 2800 at 95% CL, which 
we take as our upper limit. 

Figure 19 shows the modulation reconstructions for the 143 
GHz and 217 GHz maps separately compared to the cross esti- 
mator. The picture is complicated here by the large signals from 
noise misestimation seen in the 143 x 143 and 217 x 217 esti- 
mators, however the fact that the octopole in 143 x 143 is lower 
than in the cross-estimator indicates that the octopole signal is 
very unlikely to be mainly physical. Our measured tnl limit in 
practice represents a strong upper limit on the level of primor- 
dial tnl that could be present, since unmodelled varying small- 
scale foreground or non-constant gain/calibration would also 
only serve to increase the measured estimate compared to pri- 
mordial on average. The octopole signal does vary slightly with 
Galactic mask, though at present we cannot clearly isolate its 
origin. If more extensive analysis (for example using cross-map 
estimators at the same frequency) can identify a non-physical 
origin and remove it, the quoted upper limit on t nl would be- 
come significantly tighter. For a more extensive discussion of 
possible foreground and systematic effect issues that can af- 
fect 4-point estimators see Planck Collaboration XXVII (2013); 
Planck Collaboration XVII (2013). 

We have focussed on nearly scale invariant modulations 
here. As discussed in Planck Collaboration XXIII (2013) there 
are some interesting "anomalies" in the distribution of power in 
the Planck data, especially the hemispherical power asymmetry 
at Anax < 500. This would correspond to a r NL -like trispectrum, 
but as we have shown here the power anisotropy does not per- 
sist to smaller scales (we use £ max = 2000) except for the signal 
aligned with the CMB dipole expected from kinematic effects. 
For a primordial trispectrum to be consistent with both results 
the modulation would have to be scale-dependent on small-scale 
modes (rather than just r NL = t^(L)), so that larger small-scale 
modes are modulated more than the smallest ones. 

8. Validation of Planck results 

Here we perform a set of tests to check the robustness and sta- 
bility of our / NL measurements. As these are validation tests 
of Planck results, and not internal comparisons of bispectrum 
pipelines (already shown to be in tight agreement in Sect. 6 and 
7) we will not employ all the bispectrum estimators on each test. 
In general we choose to use two estimators on each test, in order 
to have a cross-check of the outcomes without excessive redun- 
dancy. 

8. 1 . Dependence on maximum multipole number 

The dependence on the maximum multipole number ^ max of the 
SMICA results (assuming independent shapes) is shown in Fig. 20 
(for the binned estimator) and Table 1 6 (for both the KS W and 
binned estimators). Testing the £ max dependence is easiest for the 
binned estimator, where one can simply omit the highest bins 
in the final sum when computing /nl- It is clear that we have 
reached convergence both for the values of / NL and for their error 
bars at ( max = 2500, with the possible exception of the error bars 
of the diffuse point source bispectrum. The diffuse point source 



bispectrum template is dominated by equilateral configurations 
at high L Moreover, for all the shapes except point sources, re- 
sults at ( max = 2000 are very close to those at ^ max = 2500, 
taking into account the size of the error bars. 

It is very interesting to see that at € msK ~ 500 we find a local 
/nl result in very good agreement with the WMAP-9 value of 
37.2+19.9 (Bennett et al. 2012). At these low ( max values we also 
find negative values for orthogonal /nl, although not as large or 
significant as the WMAP-9 value (which is -245 ± 100). One 
can clearly see the importance of the higher resolution of Planck 
both for the values of the different /nl parameters and for their 
error bars. 

It is also clear that the higher resolution of Planck is abso- 
lutely crucial for the ISW-lensing bispectrum; this is simply un- 
detectable at WMAP resolution. On the other hand, the high sen- 
sitivity of Planck measurements also exposes us to a larger num- 
ber of potentially spurious effects. For example we see that the 
bispectrum of point sources is also detected at high significance 
by Planck at ^ max > 2000, while remaining undetectable at lower 
resolutions. The presence of this bipectrum in the data could in 
principle contaminate our primordial /nl measurements. For this 
reason, the presence of a large point source signal has been ac- 
counted for in previous Sections by always including the Poisson 
bispectrum in a joint fit with primordial shapes. Fortunately, it 
turns out that the very low correlation between the primordial 
templates and the Poisson one makes the latter a negligible con- 
taminant for /nl, even when the residual point source amplitude 
is large. 

8.2. Dependence on mask and consistency between 
frequency channels 

To test the dependence on the mask, we have analysed the 
SMICA maps applying four different masks. Firstly the union 
mask U73 used for the final results in Sect. 7, which leaves 
73% of the sky unmasked. Secondly we used the confidence 
mask CS-SMICA89 of the SMICA technique, which leaves 89% 
of the sky. Next, a bigger mask constructed by multiplying the 
union mask U73 with the Planck Galactic mask CG60, leading 
to a mask that leaves 56% of the sky. And finally a very large 
mask, leaving only 32% of the sky, which is the union of the 
mask CL31 - used for power spectrum estimation on the raw 
frequency maps - with the union mask U73 (for mask details 
see Planck Collaboration XII 2013 for U73, CS-SMICA89, and 
CG60; Planck Collaboration XV 2013 for CL31). The results of 
this analysis are presented in Table 17 for two different esti- 
mators: binned and modal. The /nl are assumed independent 
here. In order to correctly interpret our results and conclusions, 
an important point to note is that binned results have been ob- 
tained choosing ( mix = 2500, while modal results correspond to 
f ma , = 2000. Primordial shape and ISW-lensing results and er- 
ror bars saturate at £ max = 2000 (see Sect. 8.1), so the results 
from the two estimators are directly comparable in this case. 
The Poisson (point sources) bispectrum is however dominated 
by high-f equilateral configurations and the signal for this spe- 
cific template still changes from i = 2000 to t = 2500. The 
differences in central values and uncertainties between the two 
estimators for the Poisson shape are fully consistent with the dif- 
ferent £ max values. Direct comparisons on data and simulations 
between these two estimators and the KSW estimator showed 
that Poisson bispectrum results match each other very well when 
the same ^ max is used. 

Results from the modal pipeline have uncertainties deter- 
mined from MC simulations, while the results from the binned 
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Fig. 20. Evolution of the / NL parameters (solid blue line with data points) and their uncertainties (dashed lines) for the five bispectrum 
templates as a function of the maximum multipole number £ max used in the analysis. From left to right and top to bottom the 
figures show respectively local, equilateral, orthogonal, diffuse point sources, and ISW-lensing. To better show the evolution of the 
uncertanties, they are also plotted around the final value of / NL (solid green lines without data points). The results are for SMICA, 
assume all shapes to be independent, and have been determined with the binned bispectrum estimator. 



Table 16. Results for / NL (assumed independent) of the SMICA cleaned map using different values of i 
estimators. 
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pipeline (in Table 17 and the next only) are given with Fisher 
error bars. It is very interesting to see that even with the large 
/sky = 0.32 mask, the simple inpainting technique still allows 
us to saturate the (Gaussian) Cramer-Rao bound, except for the 
ISW-lensing shape where we have a significant detection of NG 
in a squeezed configuration (so that an error estimate assuming 
Gaussianity is not good enough). Finally we note that only for 
the tests in this and in the next paragraph we adopted a faster but 
slightly less accurate version of the modal estimator than the one 
used to obtain the final /nl constraints in Sect. 7. In this faster 
implementation we use fewer modes in order to increase com- 
putational speed, and consequently we get a slight degrading of 



the level of correlation of our expanded templates with the initial 
primordial shapes. Note that the changes are small: we go from 
99% correlation for local, equilateral and orthogonal shapes in 
the most accurate (and slower) implementation to 98% correla- 
tion for equilateral and orthogonal snapes, and 95% correlation 
for local shape in the faster version. This of course still allows for 
very stringent validation tests for all the primordial shapes, and 
produces results very close to the high-accuracy pipeline, while 
at the same time increasing overall speed by almost a factor 2. 
Both versions of the modal pipeline were separately validated on 
simulations (see Sect. 6). 
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Table 17. Results for /nl (assumed independent) of the SMICA cleaned map using different masks as described in the main text 
(Sect. 8.2). Results are given for the binned and modal estimators. Uncertainties for the binned estimator in this table and the next 
are Fisher error bars. The modal estimator uses a faster and slightly less correlated expansion of the primordial templates for this 
test and the next than for other analyses (see Sect. 8.2 for more explanations). These caveats explain why the results shown in this 
table for the f s ^ y = 0.73 mask display small differences with respect to the corresponding numbers in the main results tables of 
Sect. 7, for both estimators. We also note that the binned estimator uses £ m!a = 2500 and the modal estimator f max = 2000, which 
has an impact on the point source results as explained in the main text. 



Ay = 0.89 / sky = 0.73 / sky =0.56 / sky = 0.32 

Binned 



Local 13 ±5.4 9.2 ± 5.9 11 ±6.8 6.1 ±8.9 

Equilateral 35 ± 66 -20 ± 73 -20 ± 83 39 ± 109 

Orthogonal -18 ± 36 -39 ± 39 -45 ± 45 -5 ± 59 

Diff.ps -le29 ... 14.0 ±1.3 7.7 ± 1.4 9.0 ±1.7 10.3 ± 2.2 

ISW-lensing 0.69 ± 0.26 0.91 ± 0.29 0.84 ± 0.33 0.81 ± 0.43 



Modal 

Local 12.1 ±5.5 8.4 ± 6.0 12.3 ±7.1 9.2 ± 8.7 

Equilateral 52 ± 66 -56 ± 72 -31 ±84 42 ± 104 

Orthogonal 3.3 ± 35 -31 ±40 -50 ± 47 -27 ± 59 

Diff.ps • le29 ... 20.6 ± 2.5 11.4 ±2.8 10.7 ± 3.2 12.7 ± 3.9 

ISW-lensing 0.42 ± 0.35 0.62 ± 0.40 1.1 ±0.45 0.80 ± 0.48 



Table 18. Results for / NL (assumed independent) for the raw frequency maps at 70, 100, 143, and 217 GHz with a very large mask 
(/sky = 0.32) compared to the SMICA result with the union mask U73 (/ s k y = 0.73), as determined by the binned (with ^ max = 2500) 
and modal (with £ maK = 2000) estimators. The same caveats as for the previous table (Table 17) apply here as well. 



SMICA 70 GHz 100 GHz 143 GHz 217 GHz 

Binned 



Local 9.2 ± 5.9 19.7 ± 26.0 -2.5 ± 13.2 10.4 ± 9.8 -4.7 ± 9.6 

Equilateral -20 ± 73 159 ± 188 70 ± 132 48 ±114 -9 ±114 

Orthogonal -39 ± 39 -78 ± 139 -106 ±81 -101 ±64 -84 ± 63 

Diff.ps • le29 ... 7.7 ± 1.4 (-1.4 ± 2.3)xl0 3 -4.0 ± 64 8.7 ± 6.1 14.2 ± 3.0 

ISW-lensing 0.91 ± 0.29 3.5 ± 2.2 0.35 ± 0.78 0.89 ± 0.50 0.87 ± 0.48 



Modal 

Local 8.4 ± 6.0 36.5 ± 27.2 -6.6 ± 13.6 6.6 ± 9.4 -6.5 ± 8.9 

Equilateral -56 ± 72 74 ± 193 49 ± 123 81 ± 111 28.9 ±110 

Orthogonal -31±40 -225 ± 119 -75 ± 80 -133 ±62 -112.4 ± 61 

Diff.ps -le29 11.4 ±2.8 (-2.5 ± 2.8)xl0 3 -45 ± 64 5.7 ± 7.0 25 ± 5.0 

ISW-lensing 0.62 ± 0.40 2.6 ± 2.3 0.92 ± 0.80 0.78 ± 0.60 0.85 ± 0.56 



Besides confirming again the good level of agreement be- 
tween the two estimators already discussed in Sects. 6 and 7, 
the main conclusion we draw from this analysis is that our mea- 
surements for all shapes are robust to changes in sky coverage, 
taking into account the error bars and significance levels, at least 
starting from a certain minimal mask. The / s k y = 0.89 mask is 
probably a bit too small, likely leaving foreground contamina- 
tion around the edges of the mask, though even for this mask the 
results are consistent within lcr, except for point sources (which 
might suggest the presence of residual Galactic point source con- 
tamination for the small mask). The results from the f s ^ = 0.73 
and / s ky = 0.56 masks are highly consistent. This conclusion 
does not really change when going down to f^ y = 0.32, although 
uncertainties of course start increasing significantly for this large 
mask. 

We also investigate if there is consistency between frequency 
channels when the largest mask with f^ y = 0.32 is used, and if 
these results agree with the SMICA results obtained with the com- 
mon mask. The results (assuming independent / NL ) are given 
both for the binned and the modal estimator in Table 18. As in 
the previous table, the full modal pipeline (faster but slightly less 



accurate version) has been run here, obtaining both central val- 
ues from data and MC error bars from simulations, while the 
binned pipeline (which is slower in determining full error bars 
than the modal pipeline) is used to cross-check the modal mea- 
surements and has error bars given by simple Fisher matrix es- 
timates. As one can see here and as was also checked explicitly 
in many other cases, the error bars from different estimators are 
perfectly consistent with each other and saturate the Cramer-Rao 
bound (except in the case of a significant non-Gaussian ISW- 
lensing detection). 

A detailed analysis of Table 1 8 might actually suggest that 
the agreement between the two estimators employed for this 
test, although still clearly good, is slightly degraded when com- 
pared to their performance on clean maps from different compo- 
nent separation pipelines. If we compare e.g., SMICA results in 
Table 17 to raw data results in Table 18, we see that in the for- 
mer case the discrepancy between the two estimators is at most 
of order <x/ NL /3, and smaller in most cases. In the latter case, 
however, we notice several measurements displaying differences 
of order o"/ NL /2 between the two pipelines, and the value of f^ 10 
at 70 GHz being lcr away. We explain these larger differences as 
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follows. For SMICA runs we calibrated the estimator linear terms 
using FFP6 simulations, accurately reproducing noise properties 
and correlations. On the other hand, for the present tests on raw 
frequency channels we adopted a simple noise model, based on 
generating uncorrected noise multipoles with a power spectrum 
as extracted from the half-ring null map, and remodulating the 
noise in pixel space according to the hit-count distribution. This 
approximation is expected to degrade the accuracy of the linear 
term calibration, and thus to produce a slightly lower agreement 
of different pipelines for shapes where the linear correction is 
most important. Those are the shapes that take significant con- 
tributions from squeezed triangles: local and ISW-lensing, and 
to a smaller but non-negligible extent orthogonal, i.e., exactly 
the shapes for which we find slightly larger differences. 

We conclude that no significant fluctuations are observed 
when comparing measurements from different frequency chan- 
nels (between themselves or with the clean and co-added SMICA 
map) or from different estimators on a given channel for the pri- 
mordial shapes. The same is true for the ISW-lensing shape, al- 
though it should be noted that in particular the 70 GHz chan- 
nel (like WMAP) does not have sufficient resolution to mea- 
sure either the lensing or point source contributions. The uncer- 
tainties of the point source contribution vary significantly be- 
tween frequency channels, although results remain consistent 
between channels given the error bars (when all multipoles up to 
(max = 2500 are taken into account, as performed by the binned 
estimator). This is due to the fact that this shape is dominated by 
high-^ equilateral configurations, the signal-to-noise of which 
depends crucially on the beam and noise characteristics, which 
vary from channel to channel. In the SMICA map point sources 
are partially removed by foreground cleaning, explaining the sig- 
nificantly lower value than for 217 GHz. As explained before, 
differences between the binned and modal estimators regard- 
ing point sources are due to the different values of { max (2500 
for binned and 2000 for modal), which particularly affects the 
217 GHz channel and the SMICA cleaned map. 

8.3. Null tests 

To make sure there are no hidden NG in the instrumental noise, 
we performed a set of tests on null maps. These are noise-only 
maps obtained from differences between maps with the same sky 
signal. In the first place we constructed half-ring null maps, i.e., 
maps constructed by taking the difference between the first and 
second halves of each pointing period, divided by 2. Secondly, 
we constructed a survey difference map (Survey 2 minus Survey 
1 divided by 2). A "survey" is defined as half a year of data, 
roughly the time needed to scan the full sky once; the nominal 
period of Planck data described by these papers contains two 
full surveys. Finally we constructed the detector set difference 
map ("detset 1" minus "detset 2" divided by 2). The four polar- 
ized detectors at each frequency from 100 to 353 GHz are split 
into two detector sets per frequency, in such a way that each 
set can measure all polarizations and the detectors in a set are 
aligned in the focal plane (see Planck Collaboration VI (2013) 
and Planck Collaboration XII (2013) for details on the null maps 
analysed in this Section). 

All these maps are analysed using the union mask U73 used 
for the final data results. However, in the case of the survey and 
detector set difference maps this mask needs to be increased by 
the unseen pixels. That effect only concerns a few additional pix- 
els for the detector set null map, but is particularly important for 
the survey difference map, since a survey only approximately 



covers the full sky. The final f^y of the mask used for the survey 
difference map is 64%. 

The test consisted of extracting /nl from the null maps de- 
scribed above, using only the cubic part of the bispectrum es- 
timators (i.e., no linear term correction), and keeping the same 
weights as for the full "signal + noise" analysis. This means that 
the weights were not optimized for noise-only maps, as our aim 
was not to study the bispectrum of the noise per se but rather 
to check whether the noise alone produces a three-point func- 
tion detectable by our estimators when they are run in the same 
configuration as for the actual CMB data analysis. For a similar 
reason it would have been pointless to introduce a linear term in 
this test. The purpose of the linear correction is in fact that of 
decreasing the error bars by accounting for off-diagonal covari- 
ance terms introduced by sky cuts and noise correlations when 
optimal weights are used, which is not the case here. 

Our /nl error bars for this test are obtained by running the 
estimators' cubic part on Gaussian noise simulations including 
realistic correlation properties. In the light of the above para- 
graph it is clear that such uncertainties have nothing to do with 
the actual uncertainties from CMB data, and cannot be compared 
to them. 

Since SMICA was the main component-separation method for 
our final analysis of data, we present in Table 19 the results of 
our SMICA half-ring study using the KSW, binned and modal es- 
timators, i.e., all the three main pipelines used in this paper. For 
the cleaned maps we do not have survey or detector set differ- 
ence maps. Those are, however, available for single frequency 
channels. Thus we also studied all three types of null maps for 
the raw 143 GHz channel in Table 20, using the binned estimator. 
In both tables all /nl shapes are assumed to be independent. The 
binned estimator is best suited for these specific tests as its cu- 
bic part is less sensitive to masking compared to other pipelines, 
especially modal. Therefore in this "cubic only" test, the binned 
results provide the most stringent constraints in terms of final 
error bars. 

As one can see Planck passes these null tests without any 
problems: all values found for /nl in these null maps are com- 
pletely negligible compared to the final measured results on the 
data maps, and consistent with zero within the error bars. 

8.4. Impact of foreground residuals 

In Sect. 7 we applied our bispectrum estimators to Planck 
data filtered through four different component separation meth- 
ods: SMICA, NILC, SEVEM and C-R (for a detailed descrip- 
tion of component separation techniques used for Planck see 
Planck Collaboration XII (2013)). The resulting set of /nl mea- 
surements shows very good internal consistency both between 
different estimators (as expected from our MC validation tests 
of bispectrum pipelines described in Sect. 6) and between dif- 
ferent foreground-cleaned maps. This already makes it clear that 
foreground residuals in the data are very well under control, and 
their impact on the final /nl results is only at the level of a small 
fraction of the measured error bars. In this Section we further 
investigate this issue, and validate our previous findings on data 
by running extensive tests in which we compare simulated data 
sets with and without foreground residuals from two different 
component separation pipelines, SMICA and NILC. The goal is to 
provide a MC -based assessment of the expected /nl systematic 
error from residual foreground contamination. 

For each component separation pipeline, we consider two 
sets of simulations. One set includes realistic Planck noise and 
beam, is masked and inpainted in the same way as we do for 
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Table 19. Results for /nl (assumed independent) of the SMICA half-ring null maps, determined by the KSW, binned and modal 
estimators. 



KSW Binned Modal 

SMICA half-ring 



Local -0.008 ±0.18 -0.086 ± 0.20 -0.13 ±0.35 

Equilateral -0.16 ±2.2 1.3 ±2.1 0.66 ± 2.0 

Orthogonal -0.035 ± 0.57 0.51 ± 0.57 0.14 ± 0.60 

Diff.ps • le29 ... -0.05 ± 0.60 0.03 ± 0.68 0.05 ± 0.65 

ISW-lensing (-0.06 ± 2.0)xl0~ 3 (-2.2 ± 4.7)xl0~ 3 0.009 ± 0.030 



Table 20. Results for /nl (assumed independent) of several null maps determined by the binned estimator. We consider half-ring 
(rl - r2)/2, survey (s2 - sl)/2, and detector set (dl - d2)/2 difference maps for SMICA and the raw 143 GHz channel. 



SMICA 143 GHz 143 GHz 143 GHz 

half-ring half-ring survey detector set 

Binned 



Local -0.086 ±0.20 -0.016 ± 0.073 0.43 ± 0.56 1.9 ±1.7 

Equilateral 1.3 ±2.1 3.2 ±1.8 -1.5 ± 4.2 0.9 ± 5.8 

Orthogonal 0.51 ±0.57 1.2 ± 0.6 -1.7 ±1.3 -1.3 ± 1.8 

Diff.ps • le29 ... 0.03 ± 0.68 0.19 ±1.9 3.4 ±3.2 -1.0 ±4.3 

ISW-lensing (-2.2 ± 4.7)xl0~ 3 (-0.5 ± 1.7)xlQ- 3 (-0.6 ± ll)xlQ- 3 0.033 ± 0.026 



Table 21. Summary of our /nl analysis of foreground residuals. For realistic lensed FFP6 simulations processed through the SMICA 
and NILC component separation pipelines, we report: the average /nl with and without foreground residuals added to the maps, the 
/nl standard deviation in the same two cases, and the standard deviation of the map-by-map /nl difference between the "clean" and 
"contaminated" sample. The impact of foreground residuals is clearly subdominant when compared to statistical error bars for all 
shapes. Results reported below have been obtained using the modal estimator. 



SMICA SMICA NILC NILC SMICA NILC 

clean residual clean residual residual - clean residual - clean 

Modal 



Local 7.7 ±5.9 7.8 ± 5.9 7.7 ± 5.8 7.4 ± 6.0 0.04 ±1.0 -0.27 ± 1.1 

Equilateral -0.5 ± 77 -8.7 ± 79 -0.6 ± 78 -9.0 ± 79 -8.3 ± 8.2 -8.4 ± 8.3 

Orthogonal -23 ± 41 -25 ± 41 -24 ± 40 -26 ± 41 -2.0 ± 4.7 -2.4 ± 4.8 

ISW-lensing 1.00 ±0.38 1.01 ±0.38 1.01 ± 0.38 1.02 ±0.38 0.006 ± 0.052 0.013 ± 0.052 



sequence of accidental correlations between primordial CMB 
anisotropics and foreground emission. Of course these "CMB- 
CMB-foregrounds" bispectrum terms average to zero but they 
do not cancel in the bispectrum variance 6-point function). An 
interesting point to note is that a large foreground three-point 
function will tend to produce a negative bias in the local bispec- 
trum measurements. That is because foreground emission pro- 
duces a positive skewness of the CMB temperature distribution 
("excess of photons"), and a positive skewness is in turn related 
to a negative 15 - A large negative /^J al is thus a signature 
of significant foreground contamination in the map. This is in- 
deed what we observe if we consider raw frequency maps with 
a small Galactic cut, which is why our frequency-by-frequency 
analysis in Sect. 8.2 was performed using a very large mask. For 
more details regarding effects of foreground contamination on 
primordial NG measurements in the context of the WMAP anal- 
ysis see Yadav & Wandelt 2008 and Senatore et al. 2010. 



15 While not rigorous, this argument captures the leading effect since 
Galactic foregrounds predominantly contaminate large scales. In prin- 
ciple, positively skewed, small scale foreground residuals (f > 60), or 
the negatively skewed SZ effect, can contribute positive bias. Our sim- 
ulations with foreground residuals demonstrate that these contributions 
are subdominant. 



real data, and is processed through SMICA and NILC but it 
does not contain any foreground component. The other set is 
obtained by adding to the first one a number of diffuse fore- 
ground residuals: thermal and spinning dust components; free- 
free and synchrotron emission; kinetic and thermal SZ; CO lines 
and correlated CIB. These residuals have been evaluated by ap- 
plying the component separation pipelines to accurate synthetic 
Planck datasets including foreground emission according to the 
PSM (Delabrouille et al. 2012), and are of course dependent on 
the cleaning method adopted. The simple procedure of adding 
foreground residuals to the initially clean simulations is made 
possible because we consider only linear component separation 
methods for our analysis. Linearity is in general an important 
requirement for foreground cleaning algorithms aiming at pro- 
ducing maps suitable for NG analyses. All maps in both samples 
are lensed using the LensPix algorithm. We analyse both sets 
using different bispectrum estimators (modal, KSW, binned) for 
cross-validation purposes. 

The presence of residual foreground components in the data 
can have two main effects on the measured /nl- The first is to in- 
troduce a bias in the /nl measurements due to the correlation be- 
tween the foreground and primordial 3-point function for a given 
shape. The second is to increase the error bars while leaving 
the bispectrum estimator asymptotically unbiased. This is a con- 
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Fig. 21. The measured for the first 99 maps in the lensed 
FFP6 simulation sample used for the foreground studies pre- 
sented in Sect. 8.4. We show measurements from SMICA and 
NILC processed maps both with and without foreground resid- 
uals. The horizontal solid line is the average value of the SMICA 
clean maps, and the dashed and dotted horizontal lines corre- 
spond to lcr and 2<x deviations, respectively. The plot clearly 
shows the very small impact of including residuals, and the 
very good consistency between the two component separation 
pipelines. 



In our test we built maps contaminated with foreground 
residuals by simply adding residual components to the clean 
maps. That means that the difference (f^l [dadl -f^ dn )i for the ;'-th 
realization in the simulated sample exactly quantifies the change 
in /nl due to foregrounds on that realization. In order to assess 
the level of residual foreground contamination on primordial and 
ISW-lensing shapes, first of all we consider the difference be- 
tween average values and standard deviations of /nl measured 
from the two map samples for each shape. As shown in Table 21 
we find that neither the average nor the standard deviation shows 
a significant change between the two datasets. That means that 
foreground residuals are clearly subdominant, as they do not bias 
the estimator for any shape and they do not increase the variance 
through spurious correlations with the CMB primordial signal. 

We also consider the difference / N ^ lduaI - /jS^"™ on a map-by- 
map basis and compute its standard deviation. This is used as an 
estimate of the expected bias on a single realization due to the 
presence of residuals. As already expected from the negligible 
change in the standard deviations of the two samples, the vari- 
ance of the map-by-map differences is also very small: Table 21 
again shows that it is at most about cry NL /6 for any given shape, 
where cry NL is the / NL standard deviation for that shape. As an 
example, in Fig. 21 we show the measured values of f^ 1 for 
the first 99 maps in both the SMICA and NILC samples, compar- 
ing results with and without residuals. It is evident also from this 
plot that the change in central value due to including residuals is 
very small. The very good agreement between the two compo- 
nent separation pipelines is also worth notice. 

From the study shown here and from the comparison be- 
tween different component separation methods in Sect. 7, we can 
thus conclude that the combination of foreground-cleaned maps 
and / s ky = 0.73 sky coverage we adopt in this work provide a 
very robust choice for /nl studies. 



9. Implications for early Universe physics 

The NG analyses performed in this paper show that Planck data 
are consistent with Gaussian primordial fluctuations. The stan- 
dard models of single-field slow-roll inflation have therefore 
survived the most stringent tests of Gaussianity performed to 
date. On the other hand, the NG constraints obtained on dif- 
ferent primordial bispectrum shapes (e.g., local, equilateral and 
orthogonal), after properly accounting for various contaminants, 
severely limit various classes of mechanisms for the generation 
of the primordial perturbations proposed as alternatives to the 
standard models of inflation. 

In the following Subsections, unless explicitly stated oth- 
erwise, we translate limits on NG to limits on parameters 
of the theories by constructing a posterior assuming the fol- 
lowing: the sampling distribution is Gaussian (which is sup- 
ported by Gaussian simulations); the likelihood is approximated 
by the sampling distribution but centred on the NG estimate 
(see Eisner & Wandelt 2009); uniform or Jeffreys' prior (where 
stated), over ranges which are physically meaningful, or as oth- 
erwise stated. Where two parameters are involved, the posterior 
is marginalized to give one-dimensional constraints on the pa- 
rameter of interest. 

9.1. General single-field models of inflation 

DBI models: The constraints on /n L U1 ' provide constraints on 
the sound speed with which the inflaton fluctuations can prop- 
agate during inflation. For example, for DBI models of infla- 
tion (Silverstein & Tong 2004; Alishahiha et al. 2004), where 
the inflaton field features a non-standard kinetic term, the 
predicted value of the nonlinearity parameter is fSjy = 
-(35/108)(c,T 2 - 1) (Silverstein & Tong 2004; Alishahiha et al. 
2004; Chen et al. 2007b). Although very similar to the equilat- 
eral shape, we have obtained constraints directly on the theo- 
retical (nonseparable) predicted shape (Eq. 7)). The constraint 
f™ = 1 1 ± 69 at 68% CL (see Eq. (86)) implies 

cf BI > 0.07 95% CL. (95) 

The DBI class contains two possibilities based on string con- 
structions. In ultraviolet (UV) DBI models, the inflaton field 
moves under a quadratic potential from the UV side of a warped 
background to the infrared side. It is known that this case is al- 
ready at odds with observations, if theoretical internal consis- 
tency of the model and constraints on power spectra and primor- 
dial NG are taken into account (Baumann & McAllister 2007; 
Lidsey & Huston 2007; Bean et al. 2007; Peiris et al. 2007). Our 
results strongly limit the relativistic regime of these models even 
without applying the theoretical consistency constraints. 

It is therefore interesting to look at infrared (IR) DBI mod- 
els (Chen 2005b,a) where the inflaton field moves from the IR to 
the UV side, and the inflaton potential is V((f>) = Vn - \BH 2 <p 2 , 
with a wide range 0.1 < B < 10 9 allowed in principle. In 
previous work (Beanetal. 2008) a 95% CL limit of B < 3.7 
was obtained using WMAP. In a minimal version of IR DBI, 
where stringy effects are neglected and the usual field the- 
ory computation of the primordial curvature perturbation holds, 
one finds (Chen 2005c; Chen et al. 2007b) c s =* (BN/3y\ 
n % — 1 = -4/N, where N is the number of e-folds; further, 
primordial NG of the equilateral type is generated with an 
amplitude f™ = -(35/108) [(B 2 N 2 19) - 1]. For this model, 
the range N > 60 is compatible with Planck's 3cr limits on 
« s (Planck Collaboration XXII 2013). Marginalizing over 60 < 
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N < 90, we find 



/3<0J 95% CL, 



(96) 



dramatically restricting the allowed parameter space of this 
model. 

Power-law k-inflation: These models (Armendariz-Picon et al. 
1999; Chenetal. 2007b) predict f^ n = - 170/(8 ly), where 
n s - 1 = — 3y, c 2 =! y/8. Assuming a prior of < y < 2/3, 

our constraint /^ uil = -42 ± 75 at 68% CL (see Table 9) 
leads to the limit y > 0.05 at 95% CL. On the other hand, 
Planck's constraint on n s - 1 yields the limit 0.01 < y < 
0.02 (Planck Collaboration XXII 2013). These conflicting lim- 
its severely constrain this class of models. 

Flat Models and higher derivative interactions: Flat NG can 
characterize inflationary models which arise from independent 
interaction terms different from the (n) 3 and 7r(V7r) 2 discussed 
in Sect. 2 (see also Sect. 9.2). An example is given by mod- 
els of inflation based on a Galilean symmetry (Creminelli et al. 
201 la), where one of the inflaton cubic interaction terms allowed 
by the Galilean symmetry, M 3 [7r(<9,-3 ; 7r) 2 /a 4 - IHtttt 2 + 3i/ 3 7r 3 ], 
contributes to the flat bispectrum with an amplitude = 
(35/256)(M 3 #)/(eM 2 ,) (Creminelli et al. 2011a). Here, n is the 
relevant inflaton scalar degree of freedom, e the usual slow-roll 
parameter and a the scale factor and H the Hubble parame- 
ter during inflation. Our constraint = 37 ± 77 at 68% CL 
(see Table 11) leads to (M 3 //)/(eM 2 ,) = 270 ± 563 at 68% 
CL. These interaction terms are similar to those arising in gen- 
eral inflaton field models that include extrinsic curvature terms, 
e.g., parameterized in the Effective Field Theory approach as 
M 2 A(dij7r) 2 /a 4 (Bartolo et al. 2010a), which contribute to a flat 
bispectrum with an amplitude = (50/108) M 2 /(c 2 eMpj). In 
this case, we obtain M 2 /(c 2 eM 2 ,) = 80 ± 166 at 68% CL. 



9.2. Implications for Effective Field Theory of Inflation 

The effective field theory approach to inflation (Weinberg 2008; 
Cheung et al. 2008) provides a general way to scan the NG pa- 
rameter space of inflationary perturbations. For example, one 
can expand the Lagrangian of the dynamically relevant degrees 
of freedom into the dominant operators satisfying some under- 
lying symmetries. We will focus on general single-field models 
parametrized by the following operators (up to cubic order) 
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Fig. 22. 68%, 95%, and 99.7% confidence regions in the param- 
eter space (f^, 



./SJ* ), defined by thresholding^- 2 as described 



in the text. 
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where n is the scalar degree of freedom (£ = -Hji). The mea- 
surements on /nl"' 1 an d /ni* can ^ e usec ' t0 constram me mag- 
nitude of the inflaton interaction terms 7r(<9,7r) 2 and (n) 3 which 
give respectively f^ 11 = -(85/324)(c.- 2 - 1) and f^ 2 = 

-(10/243)(c" 2 - 1) [c 3 + (3/2)c 2 ] (Senatore et al. 2010, see also 
Chen et al. 2007b; Chen 2010b). These two operators give rise 
to shapes that peak in the equilateral configuration that are, 
nevertheless, slightly different, so that the total NG signal will 
be a linear combination of the two, possibly leading also to 
an orthogonal shape. There are two relevant NG parameters, 
c s , the sound speed of the the inflaton fluctuations, and M 3 
which characterizes the amplitude of the other operator 7f 3 . 



Fig. 23. 68%, 95%, and 99.7% confidence regions in the single- 
field inflation parameter space (c s , c 3 ), obtained from Fig. 22 via 
the change of variables in Eq. (98). 



Following Senatore et al. (2010) we will focus on the dimension- 
less parameter c 3 (Cj7 2 - 1) = 2M 4 cll{HM 2 ^}. For example, DBI 
inflationary models corresponds to c 3 = 3(1 - c 2 )/2, while the 
non-interacting model (vanishing NG) correspond to c s = 1 and 
M 3 =0 (or c 3 (c- 2 - 1) = 0). 

The mean values of the estimators for equilateral and orthog- 
onal NG amplitudes are given in terms of c s and c 3 by 



/.equil 



/-ortho 
/NL 



1 -c 2 

-(-0.275 + 0.0780A) 



1-c 2 



-(0.0159 - 0.0167A) 



(98) 
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where A = — (c* + (2/3)c^), and the coefficients are com- 
puted from the Fisher correlation matrix between the equilat- 
eral and orthogonal template bispectra and the theoretical bis- 
pectra arising from the two operators 7r(V7r) 2 and A 3 . Given 
our constraints on /S ml and f^ 10 , and the covariance matrix 
C of the joint estimators, we can define a x 2 statistic given by 
X 2 (c3,c s ) = « r (c3,c s )C _1 u(c3,c s ), where the vector u is given 
by y'(c3,c s ) = f'(cj,,c s ) - f' p . f' p , where /^{equilateral, orthogo- 
nal), are the joint estimates of the equilateral and orthogonal /nl 
measured by Planck and /'(C3, c s ) is given by Eq. (98). Figure 22 
shows the 68%, 95%, anc/99.7% confidence regions for /S^ and 
/°£ ho , obtained by requiring^ 2 < 2.28,5.99, and 11.62 respec- 
tively, as appropriate for a^ 2 variable with two degrees of free- 
dom. The corresponding confidence regions in the (c3,c s ) pa- 
rameter space are shown in Fig. 23. After marginalizing over C3 
we find the following conservative bound on the infiaton sound- 
speed 

c s > 0.02 95% CL . (99) 

Note that we have also looked explicitly for the non-separable 
shapes in Sect. 7.3.1, in particular the two effective field theory 
shapes and the DBI inflation shape (see Eqs. (5, 6, 7)) . 
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9.3. Multi-field models 

Curvaton models: Planck NG constraints have interesting im- 
plications for the simplest adiabatic curvaton models. They pre- 
dict (Bartolo et al. 2004d,c) 
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(100) 



for a quadratic potential of the curvaton field (Lyth & Wands 
2002; Lythetal. 2003; Lyth & Rodriguez 2005; Sasaki et al. 
2006), where r D = [3p curva ton/(3pcui-vaton + 4p radiation )] D is the 
"curvaton decay fraction" evaluated at the epoch of the curva- 
ton decay in the sudden decay approximation. Assuming a prior 
< r D < 1, given our constraint /^° L cal = 2.7 + 5.8 at 68% CL, 
we obtain 

r D > 0.15 95% CL. (101) 

In Planck Collaboration XXII (2013) a limit on is derived 
from the constraints on isocurvature perturbations under the as- 
sumption that there is some residual isocurvature fluctuations in 
the curvaton field. For this restricted case, they find r^ > 0.98 
(95% CL), compatible with the constraint obtained here. 

Quasi-single field inflation: It is beyond the scope of this pa- 
per to perform a general multi-field analysis employing the local 
NG constraints. However, we have performed a detailed anal- 
ysis for the quasi-single field models (see Eq. (12)). Quasi- 
single field (QSF) inflation models (Chen & Wang 2010b,a; 
Baumann & Green 2012) are a natural consequence of inflation 
model-building in string theory and supergravity (see Sect. 2.2). 
In addition to the infiaton field, these models have extra fields 
with masses of order the Hubble parameter, which are stabilized 
by supersymmetry. A distinctive observational signature of these 
massive fields is a one-parameter family of large NG whose 
squeezed limits interpolate between the local and the equilat- 
eral shape. Therefore, by measuring the precise momentum- 
dependence of the squeezed configurations in the NG, in prin- 
ciple, we are directly measuring the parameters of the theory 
naturally determined by the fundamental principle of supersym- 
metry. These models produce a bispectrum (Eq. (12)) depending 
osi 

on two parameters v, ni , with a shape that interpolates between 



Fig. 24. 68%, 95%, and 99.7% confidence intervals for v and 
osi 

/nl for quasi-single field inflation. The best fit value of v = 1.5, 
/nl 1 = 4.75 is marked with an X. The contours were calculated 
using MC methods by creating 2 x 10 9 simulations using the ft 
covariance matrix around this best fit model. 



the local shape, where v = 1.5 and the equilateral shape, where 
v = 0. 

Results are shown in Fig. 24 (see Sect. 7.3.6 for details of 
the analyses). The best fit value corresponds to v = 1.5, /nl = 
4.79 which would imply, within this scenario, that the extra field 
different from the infiaton has a mass m «: H. However, the 
figure shows that there is no preferred value for v with all values 
allowed at 3<x. 

Alternatives to inflation: Perhaps the most striking example is 
given by the ekpyrotic/cyclic models (for a review, see Lehners 
2010) proposed as alternative to inflationary models. Typically 
they predict a local NG |/^ al | > 10. In particular, the so-called 
"ekpyrotic conversion" mechanism (in which isocurvature per- 
turbations are converted into curvature perturbations during the 
ekpyrotic phase) yields f^ dX = -(5/12) c 2 , where c\ is a param- 
eter in the potential, requiring 10 > cj > 20 for compatibility 
with power spectrum constraints. This case was ~ 4<x discrepant 
with WMAP data, and with Planck it is decisively ruled out given 
our bounds /^° L cal = 2.7 + 5.8 at 68% CL (see Table 9) yield- 
ing c\ < 4.2 at 95% CL. The predictions for the local bispec- 
trum from other ekpyrotic models (based on the so called "ki- 
netic conversion" taking place after the ekpyrotic phase) yield 



f^f = (3/2) * 3 Ve+5, where e ~ 100 is natural (Lehners 2010). 
Assuming a prior -1 < K3 < 5, we obtain -0.8 < K3 < 0.5 at 
95% CL, dramatically restricting the viable parameter space of 
this model. 



9.4. Non-standard inflation models 

Constraints on excited initial states: Results from Sect. 7.3 con- 
strain a variety of models with flattened bispectra, often in 
combination with a non-trivial squeezed limit. The most no- 



45 



Planck Collaboration: Planck 2013 Results. XXIV. Constraints on primordial NG 



table examples are bispectra produced in excited initial state 
models (non-Bunch-Davies vacua), which can be generated by 
strong disturbances away from background slow-roll evolu- 
tion or additional trans-Planckian physics (Chen et al. 2007b; 
Holman & Tolley 2008; Meerburg et al. 2009; Agullo & Parker 
2011). The constraints we have obtained are summarized in 
Table 1 1, and cover four representative cases (see Eq. (14, 15)) 
in the literature. We find no strong evidence for these flattened 
bispectra in the Planck data after subtraction of the ISW-lensing 
signal, with which all these models have some correlation. This 
is consistent with an earlier constraint on the NBD model from 
WMAP (Fergusson et al. 2012). However, this investigation is 
limited by the present resolution of the polynomial modal es- 
timator (n max = 600), so more strongly flattened models are not 
excluded by this analysis. 

Directional dependence motivated by fields: Directionally- 
dependent bispectra (Eq. (19)), motivated by inflation with 
gauge fields, have also been constrained (see Table 11). For ex- 
ample, models with a kinetic term of the vector field(s) £, = 
-I 2 (<p)F 2 /4, where F 2 is the strength of the gauge field, and /(</>) 
is a function of the inflaton field which, with an appropriate time 
dependence (see, e.g., Ratra 1992), can generate vector fields 
during inflation. For these models the L — 0, 2 modes in the 
bispectrum are excited with f^ L = Xi(\g m \/0A) (Nk 3 /60), where 
X L=0 = (80/3) and X L=2 = -(10/6), respectively (Barnaby et al. 
2012a; Bartolo et al. 2013; Shiraishi et al. 2013). Here g t is the 
amplitude of a quadrupolar anisotropy in the power spectrum 
(see, e.g., Ackerman et al. (2007)) and N is the number of e- 
folds before the end of inflation when the relevant scales exit 
the horizon. These modes therefore relate the bispectrum am- 
plitude to the level of statistical anisotropy in the power spec- 
trum. Marginalizing over a prior 50 < TV < 70 assuming uni- 
form priors on g*, our constraints in Table 1 1 lead to the limits 
-0.05 <g t < 0.05 and -0.36 < g* < 0.36 from the L = 0, L = 2 
modes respectively (95% CL). Note, however, that in the current 
analysis only a modest correlation was possible with the shape 
corresponding to the L — 2 mode. These results apply to all mod- 
els where curvature perturbations are sourced by a I 2 ((p)F 2 term 
(see references in Shiraishi et al. 2013). 

Feature models: Non-scale-invariant oscillatory bispectrum 
shapes can be generated by sharp or periodic features in the 
inflaton potential, with particular recent interest on axion mon- 
odromy models motivated by string theory (see Sect. 2.3). We 
have undertaken a survey of simple feature models (Eq. (16)) 
which have a periodicity accessible to the polynomial modal 
estimator (wavenumbers K = k\ + k 2 + k^ > 0.01), a two- 
parameter space spanned by K and a phase </>. There are inter- 
esting best fit models for the Planck CMB bispectrum around 
K = 0.01875, = (that is, with a large-^ bispectrum pe- 
riodicity around At = 260), with results shown in Table 12. 
We note important caveats on the statistics of parameter sur- 
veys like this in Sect. 7.3.3; given the large numbers of feature 
models studied, we have to apply a higher threshold for statis- 
tical significance as shown for a survey of 200 Gaussian sim- 
ulations. This feature survey takes forward earlier results for 
the WMAP data (Fergusson et al. 2012), with the apparent fit 
reflecting the signal observed in the Planck CMB reconstruc- 
tion (see Fig. 7). Most attention on feature models has been 
motivated by the simplest single-field case for which there are 
correlated signals predicted in the bispectrum and power spec- 
trum (see e.g., Chen et al. 2007a; Adshead et al. 2012). In this 
case, regions with small k ~ 0.001 are favoured for producing 



an observable bispectrum , given existing WMAP power spec- 
trum constraints on these models. Periodicities A{ < 20 are 
anticipated (see Adshead et al. 2012) which are not accessible 
to the present modal bispectrum estimator analysis, but which 
are discussed in (Planck Collaboration XXII 2013). Conversely, 
the Planck feature model survey using the power spectrum 
(Planck Collaboration XXII 2013) does cover bispectrum scales 
and parameters investigated in this paper. An analysis of 
the Bayesian posterior probability (Planck Collaboration XXII 
2013) does not appear to provide evidence favoring parameters 
associated with the current best-fit bispectrum model. More de- 
tailed analysis using the specific bispectrum envelope for the 
single-field feature solution is being undertaken. 

Resonance models: We have also investigated resonance models 
of Eq. (17) such as axion monodromy and enfolded resonance 
models of Eq. (18), in which non-Bunch-Davies vacua are ex- 
cited by sharp features. This limited analysis focuses on period- 
icities associated with the best-fit feature model and the results 
are described in Tables B.l and B.2. No significant signal was 
found in this domain for either of these two models. However, 
we note that the logarithmic dependence of the bispectrum cre- 
ates challenges at low k, as we must ensure important features 
do not fall below the modal resolution. This restricts the present 
survey range, which will be extended in future. Again, we note 
that most attention on these models has focused on higher t- 
periodicities than those accessible to the present modal estima- 
tor (see e.g., Flauger & Pajer 2011; Peiris et al. 2013), for the 
same reason as the feature models. A resonance model survey 
using the Planck power spectrum has been undertaken and the 
results can be found in Planck Collaboration XXII (2013); how- 
ever, this also currently excludes high frequencies that have re- 
ceived attention in the literature. 

Warm inflation: This model, where dissipative effects are impor- 
tant, predicts = -151n(l + r d /14) - 5/2 (Moss & Xiong 
2007) where the dissipation parameter = r/(3H) must be 
large for strong dissipation. The limit from WMAP is < 2.8 x 
10 4 (Moss & Xiong 2007). Assuming a prior < log 10 r d < 4, 
our constraint / ] ^ rmS = 4 ± 33 at 68% CL (see Sect. 7.3.5) 
yields a limit on the dissipation parameter of log 10 rd < 2.6 
(95% CL), improving the previous limit by nearly two orders 
of magnitude. The strongly-dissipative regime with r& > 2.5 still 
remains viable; however, the Planck constraint puts the model 
in a regime which can lead to an overproduction of gravitinos 
(see Hall & Peiris (2008) and references within). 

9.5. Implications of CMB thspectrum results 

The non-detection of local-type / NL by Planck raises the imme- 
diate question of whether there might still be a large trispec- 
trum. In this first analysis we have focused on the local shape 
t nl- Most inflation models predict t nl ~ 0(f 2 L ) (Byrnes et al. 
2006; Elliston et al. 2012), and hence given our tight /nl limits, 
would be predicted to be very small. This is easily consistent 
with our conservative upper limit t nl < 2800, and also with 
the significantly smaller signals found in the modulation dipole 
and quadrupole. Our upper limit also restricts the freedom of 
curvaton-like models with quadratic terms that are nearly uncor- 



16 Planck Collaboration XVI (2013) confirms an anomaly in the 
power spectrum at 20 < I < 40 first noted in WMAP, which leads to an 
improvement in likelihood when fitted with a feature in the inflationary 
potential (Peiris et al. 2003). Unfortunately, the best-fit parameters in 
this case do not lead to an observable bispectrum (Adshead et al. 201 1). 
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related with the curvature perturbation, which could in principle 
have /nl ~ but a large trispectrum (Byrnes & Choi 2010). 

10. Conclusions 

In this paper we have derived constraints on primordial NG, us- 
ing the CMB maps derived from the Planck nominal mission 
data. Using three optimal bispectrum estimators - KSW, binned, 
and modal - we obtained consistent values for the primordial 
local, equilateral, and orthogonal bispectrum amplitudes, quot- 
ing as our final result /^° L cal = 2.7 + 5.8, f^ il = -42 + 75, 

and f° T t ho = -25 ± 39 (68 % CL statistical). Hence there is 
no evidence for primordial NG of one of these shapes. We did, 
however, measure the ISW-lensing bispectrum expected in the 
ACDM scenario, as well as a contribution from diffuse point 
sources, and these contributions are clearly seen in the form 
of the associated skew-Q. These results have been confirmed 
by measurements using the wavelet bispectrum, and Minkowski 
functional estimators, and demonstrated to be stable for the four 
different component separation techniques SMICA, NILC, SEVEM, 
and C-R, showing their robustness against foreground residuals. 
They have also passed an extensive suite of tests studying the 
dependence on the maximum multipole number and the mask, 
consistency checks between frequency channels, and several null 
tests. In addition, we have summarized in this paper an extensive 
validation campaign for the three optimal bispectrum estimators 
on Gaussian and non-Gaussian simulations. 

Extending our analysis beyond estimates of individual 
shape amplitudes, we presented model-independent, three- 
dimensional reconstructions of the Planck CMB bispectrum us- 
ing the modal and binned bispectrum estimators. These results 
were also shown to be fully consistent between the different 
component separation techniques even for the full bispectrum, 
and contained no significant NG signals of a type not captured 
by our standard templates at low multipole values. At high mul- 
tipoles, some indications of unidentified NG signals were found, 
as also evidenced by the results from the skew-Q estimator. 
Further study will be required to ascertain whether these indi- 
cations are due to foreground residuals, beams, data processing, 
or a more interesting signal. 

Using the modal decomposition, we have presented con- 
straints on key primordial NG scenarios, including general 
single-field models with non-separable shapes, excited initial 
states (non-Bunch-Davies vacua), and directionally-dependent 
vector models. We have also undertaken an initial survey of 
scale-dependent feature and resonance models. 

Moving beyond three-point correlations, we have obtained 
limits from the Planck data on the amplitude of the local four- 
point function. Our trispectrum reconstruction yielded a signal 
consistent with zero except for an anomalously large octopole. 
Frequency dependence indicated that this was unlikely to be pri- 
mordial, but allowing the signal to be primordial we placed an 
upper limit r NL < 2800 (95% CL). 

We have discussed the impact of these results on the infla- 
tionary model space, and derived bispectrum constraints on a se- 
lection of specific inflationary mechanisms, including both gen- 
eral single-field inflationary models (extensions to the standard 
single-field slow-roll case) and multi-field models. Our results 
led to a lower bound on the speed of sound, c s > 0.02 (95% 
CL), in an effective field theory approach to inflationary models 
which includes models with non-standard kinetic terms. Strong 
constraints on other scenarios such as IR DBI, ^-inflation, infla- 
tionary models involving gauge fields, and warm inflation have 



been obtained. Within the class of multi-field models, our mea- 
surements limit the curvaton decay fraction to ru > 0.15 (95% 
CL). Ekpyrotic/cyclic scenarios were shown to be under severe 
pressure from the Planck data: we robustly ruled out the so- 
called "ekpyrotic conversion" mechanism, and found that the pa- 
rameter space of the "kinetic conversion" mechanism is severely 
limited. 

With these results, the paradigm of standard single-field 
slow-roll inflation has survived its most stringent tests to-date. 
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Appendix A: Expected level of agreement from 
bispectrum estimators with correlated weights 

The estimator cross-validation work presented in Sect. 6.1 was 
based on comparing results from different estimators using sets 
(typically 50 to 100 simulations in size) of both Gaussian and 
non-Gaussian simulations. We started from idealized maps (e.g., 
full-sky, noiseless, Gaussian simulations) and then went on to 
include an increasing number of realistic features at each addi- 
tional validation step. This allowed better testing and character- 
ization of the response of different pipelines and bispectrum de- 
compositions to various potential spurious effects in the data. As 
a preliminary step, we derived a general formula describing the 
expected level of agreement between estimators with different 
but strongly correlated weights, with the simplifying assumption 
of full-sky measurements and homogeneous noise. This theoret- 
ical expectation, summarized by Eq. (84), was then used as a 
benchmark against which to assess the quality of the results. The 
aim of this appendix is to describe in detail how we obtained 
Eq. (84). 

Let us consider the idealized case of full-sky, noiseless CMB 
measurements (note that the following conclusions also work 
for homogeneous noise, because the pure cubic /nl estimators 
without linear term corrections are still optimal) . Under these 
assumptions, the / NL estimator for a given CMB shape B^ 2 f 3 
can be written simply as (see Sect. 3 for details): 



1 

F 



I E 



M^2^3^ ]mi ^2^2 ^3 m 3 



(A.l) 



where B* £ ( is the angle-averaged bispectrum template for a 
given theoretical shape, ae imi ae 2m2 at 2m is a bispectrum estimate 
constructed from the data, and F is the normalization of the es- 
timator, provided by the Fisher matrix number 
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As explained in Sect. 3, the different /nl estimation techniques 
implemented in this work can be seen as separate implementa- 
tions of the optimal estimator of Eq. (A.l). Each implementa- 
tion is based on expanding the angular bispectrum as a sum of 
basis templates defined in different domains: for example in our 
analyses we built templates out of products of wavelets at dif- 
ferent scales, cubic combinations of 1 -dimensional polynomials 
and plane waves (what we call "modal estimator" in the main 
text), and ^-binning of the bispectrum (what we called "binned 
estimator" in the main text). Our initial theoretical templates 



50 



Planck Collaboration: Planck 2013 Results. XXIV. Constraints on primordial NG 



in this work are the local, equilateral and orthogonal separable 
cases used in the KSW and skew-Q estimators. In this sense 
the KSW/skew-CV estimator provides an "exact" estimate of /nl 
for this choice of shapes, while the other pipelines provide an 
approximate result that approaches KSW measurements as the 
expansions get more accurate. The differences between results 
from different pipelines became smaller and smaller as the ap- 
proximate expanded templates converge to the starting one (e.g., 
by increasing the number of ^-bins or the number of wavelets 
and polynomial triplets). The degree of convergence can be mea- 
sured through the correlation coefficient r between the initial bis- 
pectrum and its reconstructed expanded version. The correlation 
coefficient is, as usual, defined as 



y L r2 L 3 L r2'3 



(A.3) 
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where the label "th" denotes the initial bispectrum shape to fit to 
the data, and "exp" is the approximate expanded one. The corre- 
lator between shapes naturally defines a scalar product in bispec- 
trum space, and in the following the operation of correlating two 
shapes will be denoted by the symbol ( , ), so that the definition 

above would read r = <B th , fi exp >/ V<B exp , B exp )<B th , B th ). 

Whichever basis and separation scheme we have chosen, let 
us call Uniti, €2, (3) the n-th template in the adopted bispectrum 
expansion, and a„ the coefficients of the expansion, so that we 
can write: 



(A.4) 



n=0 



From now on we will also call ^,,(^1,^2,^3) the "modes" of our 
expansion, N exp is the number of modes we are using to approxi- 
mate the starting template in order to obtain a correlation coeffi- 
cient r. We will also assume that the modes form an orthonormal 
basis, that is: 

(K m ,<R n - l )=$Z< ( A -5) 

where §"{ x is the Kronecker delta symbol. The orthogonality as- 
sumption does not imply loss of generality since any starting 
set of modes can always be rotated and orthogonalized. We now 
consider an expansion with a number of coefficients A^h > N exp 
that perfectly reproduces the initial bispectrum (i.e., r = 1), and 
is characterized by the same modes and coefficients as the pre- 
vious one up to iVexp- So we can write 



Mi, 



(A.6) 



We now build two optimal estimators of /nl for the shape B±: 
an "exact" estimator and an "approximate" one. In the exact es- 
timator we fit the actual B± shape to the data and obtain the esti- 
mate / t h, while in the approximate estimator we fit the expanded 
shape Z? exp to obtain / exp . We want to understand by how much 
the "exact" and "approximate" / NL measurements are expected 
to differ, as a function of the correlation coefficient r between the 
weights B t h and B exp that appear in the two estimators. 

For each mode template % n {i\ , t% t (3) we can build an opti- 
mal estimator (following Eq. (A.l)) in order to fit the mode to 
the data and get a corresponding amplitude estimate B„. In other 
words, the observed bispectrum can then be reconstructed as in 
Eq. (A.4), but with coefficients B„ instead of a„. Due to the or- 
thonormality of the %„ the B coefficients have unit variance. It 



is then possible to show (Fergusson et al. 2010a) that the /nl es- 
timate for a given B t h with expansion parameters a n is given by 



(A.l) 



In the light of all the above, the exact and approximate estimators 
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Thanks to the orthonormality properties of the modes, we can 
easily relate the Fisher matrix normalization F to the expansion 
coefficients a: 
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In analogous fashion we can derive an expression for the correla- 
tion coefficient r between the two estimators we are comparing. 
It is easy to check that the following holds: 



y«exp 2 



(All) 



Y N,h a 2 

and using the equation just derived above we can also write r 2 = 
F exp /F±, i.e., the square of the correlation coefficients between 
the estimators equals the ratio of the normalizations. 

Armed with this preliminary notation we can now calcu- 
late the expected scatter between the exact and approximate /nl 
measurement when we apply the two estimators to the same 
set of data. In order to quantify it we will consider the variable 
<5f NL = f^ L - /^ xp , and calculate its standard deviation. We find 

,2 
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(A. 12) 

where we made use of Eq. (A.l 1). It is easy to show that or- 
thonormality of the % templates implies no correlation of the 
amplitudes B, i.e., {B p B q = 5 q p ). A straightforward calculation 
then yields: 



= ^2 a "- 2r2 2 a » +r4 2^- (A. 13) 
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This, together with Eqs. (A. 10, A.l 1) finally gives: 

Vl -r 2 

o- gfNL = Ath , (A. 14) 

where A t h is the standard deviation of the exact estimator. Eq. 
(A. 14) is the starting point of our validation tests. It provides 
an estimate of the expected scatter between /nl estimators with 
correlated weights, as a function of the /nl error bars and of 
the correlation coefficient r. Note that this formula has been ob- 
tained under the simplifying assumptions of Gaussianity, full- 
sky coverage and homogeneous noise. For applications dealing 
with more realistic cases we might expect the scatter to become 
larger than this expectation, while remaining qualitatively con- 
sistent. 

In our tests we started from sets of simulated maps and com- 
pared /nl results from different pipelines and shapes on a map- 
by-map basis. In our validation tests the correlation levels be- 
tween templates in different expansions were varying between 
r ~ 0.95 and r ~ 0.99, depending on the estimators and the 
shapes under study. Using the formula above, this corresponds 
to an expected scatter 0.15 A < <x < 0.3 A, where A is the /nl 
standard deviation for the shape under study. 

In Sect. 6.1 we presented several applications to simulated 
data sets, showing that the recovered results are fully consistent 
with these expectations, and thus the different pipelines are fully 
consistent with each other. 

Appendix B: Targeted Bispectrum Constraints 

This Appendix provides further tables of results for primordial 
models, extending those given in Sect. 7.3, notably for reso- 
nance models, while it also gives additional validation checks for 
the modal bispectrum estimator, beyond the extensive tests de- 
scribed already in the paper for the standard bispectrum shapes. 
In particular, using each of the SMICA NILC and SEVEM maps, 
we will quote results for the main paradigms for non-standard 
bispectrum models, including comparisons for the feature model 
results. Remarkably consistent results are again obtained using 
the different foreground-separation methodologies and using dif- 
ferent modal basis functions. 

B. 1 . Additional results for resonance models and enfolded 
resonance models 

As described in Sect. 7.3, we have surveyed resonance mod- 
els (Eq. (17)) in the region of most interest for feature models, 
that is, with comparable periodicities to those with described in 
Table 12 near the best-fit feature model. The results from this 
initial survey for the SMICA component-separated map produced 
no significant signal, with the results Table B.l. We note that 
while the feature and resonance models proved similar for the 
WMAP analysis (Fergusson et al. 2012), wavelengths are much 
more differentiated for Planck and so it can be difficult to re- 
solve the shortest wavelengths at low £. This is the key limitation 
on surveys with the present estimator and will be circumvented 
in future, by using specific separable templates to improve the 
overall resolution. Just as local and ISW templates can be incor- 
porated into the analysis, so can targeted feature modes. Note 
that consistent results were obtained using the NILC and SEVEM 
maps, though we only show results here for feature models with 
an envelope (see discussion below with Table B.4). 

The enfolded resonant (or non-Bunch-Davies feature) mod- 
els (Eq. (18)) were also surveyed for these periodicities and also 



yielded no significant signal; see Table B.2. These are interest- 
ing shapes which hold out the prospect of exhibiting both oscil- 
latory and flattened features observed in the Planck bispectrum 
reconstruction, see Fig. 7. Due to similar resolution restrictions, 
only relatively large multipole periodicities (t > 100) have been 
surveyed in the present work, again searching around the peri- 
odicities exhibited for feature models. 



B.2. Comparison of f Nh results from SMICA, NILC and SEVEN 
foreground-cleaned maps 

Tables B.3 and B.4 compare bispectrum results extracted from 
Planck maps created using the three different component- 
separation techniques, SMICA, NILC and SEVEM. In addition 
to the models discussed in Sect. 7.3, we have also included 
the constant model which is defined by B cm!il {k\,k2,k's) = 
6A 2 f™ st /(kik 2 k 3 ) 2 . The abbreviations EFT denotes the effective 
field theory single field shapes, NBD the non-Bunch-Davies (ex- 
cited initial state) models, vector models are gauge field inflation 
with directional dependence, along with DBI, Ghost and Warm 
inflation, also described previously. The expression "cleaned" 
refers to removal of the predicted ISW-lensing signal and the 
measured point source signal, unless stated otherwise. 

We note that there is good consistency between the differ- 
ent foreground-separation techniques for all models, whether 
equilateral, flattened, or squeezed as shown in Table B.3. For 
the non-scaling case, differences for the best-fit feature models 
were below l/3cr confirming the interesting results discussed in 
Sect. 3.2.3, see Table B.4. 

B.3. Comparison of / N l results from ISW Fourier basis with 
hybrid polynomials 

As described in Sect. 3.2.3, the modal bispectrum estima- 
tor can flexibly incorporate any suitable basis functions with 
which to expand the bispectrum and separably filter CMB 
maps (Fergusson et al. 2010a). For the Planck analysis, we have 
evolved two sets of basis functions to fulfil three basic crite- 
ria: first, to provide a complete basis for the bispectrum up to 
a given ^-resolution, secondly, to accurately represent primor- 
dial models of interest and, thirdly, to incorporate the CMB 
ISW-lensing signal, which with diffuse point sources provides 
a significant secondary signal which must be subtracted. The 
first basis functions are « max = 600 polynomials (closely re- 
lated to Legendre polynomials) which are supplemented with the 
Sachs- Wolfe local mode in order to more accurately represent 
the squeezed limit; enhanced orthogonality is preconditioned by 
choosing these from a larger set of polynomials. The second ba- 
sis functions are n max = 300 Fourier modes, augmented with 
the same SW local mode, together with the the separable ISW- 
lensing modes. Both modal expansions proved useful, provid- 
ing important validation and cross-checks, however, the twofold 
resolution improvement from the polynomials meant that most 
quantitative results employed this basis. This improved resolu- 
tion was particularly important in probing flattened models on 
the edge of the tetrapyd. 

In Fig. 7, we show a direct comparison between the modal 
reconstructtion of the 3D bispectrum using the polynomial and 
Fourier mode expansions. The basic features of the two recon- 
structions are in good agreement, confirming a central feature 
which changes sign at low t and a flattened signal beyond as 
discussed previously in Sect. 7.2. Notably the polynomial ba- 
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Table B.l. Results from a limited /nl survey of resonance models of Eq. (17) with 0.25 < k c < 0.5 using the SMICA component- 
separated map. These models have a large-^ periodicity similar to the feature models in Table 12. 



Phase = = tt/5 = 2n/5 = 3n/5 = 4^/5 = n 

Wavenumber / NL ± A/ NL / NL ± A/ NL / NL ± A/ NL / NL ± A/ NL / NL ± A/ NL / NL ± A/ NL 



fc c = 0.25 -16 ±57 6 ±63 19 ±67 31 ± 69 38 ± 68 -6 ± 60 

fc c =0.30 -66 ±73 -57 ±74 -44 ± 73 -26 ± 72 -7 ±71 -65 ± 73 

fc c =0.40 5 ±57 40 ±66 55 ±71 63 ± 73 63 ±71 22 ±61 

k c = 0.45 25 ± 56 34 ± 59 36 ± 62 34 ± 67 27 ± 69 30 ± 56 

fc c = 0.50 -2 ±65 -13 ±72 -16 ±69 -16 ±60 -14 ±55 -7 ±71 



Table B.2. Results from a limited /nl survey of non-Bunch-Davies feature models (or enfolded resonance models) of Eq. (18) with 
4 < k c < 12, again overlapping in periodicity with the feature model survey. 



Phase = <f> = n/4 = ^/2 = 3^/4 

Wavenumber /nl ± A/ NL (cr) / NL ± A/ NL (cr) / NL ± A/ NL (cr) / NL ± A/ NL (cr) 



k c = 4 11 ± 146 ( 0.1) 2±145( 0.0) -7 ± 143 ( 0.0) -15 ± 142 (-0.1) 

kc = 6 52 ± 202 ( 0.3) 63 ± 203 ( 0.3) 72 ± 201 ( 0.4) 80 ± 197 ( 0.4) 

k c = 8 100 ± 190 ( 0.5) 130 ± 189 ( 0.7) 158 ± 189 ( 0.8) 183 ± 190 ( 1.0) 

k c = 10 188 ± 241 ( 0.8) 210 ± 242 ( 0.9) 230 ± 242 ( 1.0) 248 ± 243 ( 1.0) 

k c = 12 180 ± 307 ( 0.6) 171 ± 310 ( 0.6) 158 ± 312 ( 0.5) 142 ± 314 ( 0.5) 



sis, with double the resolution, preserves the large-scale features 
observed in the Fourier basis. 

In Table B.5, results from both basis expansions are shown 
for a variety of non-separable models. These demonstrate con- 
sistent results where the Fourier basis had sufficient resolution, 
as indicated by the ratio of the variance reflecting the Fisher 
ratio (i.e., above 90% correlation as indicated by the results in 
Appendix A). Note that we independently determine the estima- 
tor correlation between the exact solution and primordial decom- 
position and then at late times with the CMB decomposition; we 
use a polynomial basis as the overall benchmark here. This anal- 
ysis also includes several feature models (phase = 0) show- 
ing good agreement from the Fourier basis while the Fisher esti- 
mates remain accurate. Again, the hybrid Fourier basis degrades 
in accuracy towards k = 0.02 as it reaches its resolution limit, 
when the variance disparity rises towards 10%. With n max = 600 
modes and { max = 2000, the polynomial basis retains a good 
correlation for all primordial feature models for k > 0.01. The 
accuracy and robustness of the feature model results have been 
verified using f max = 1500 for the polynomial expansion, for ex- 
ample, obtaining 3. lcr with the SMICA map for the best-fit model 
(K = 0.01875, = 0). 



APC, AstroParticule et Cosmologie, Universite Paris Diderot, 
CNRS/IN2P3, CEA/lrfu, Observatoire de Paris, Sorbonne Paris 
Cite, 10, rue Alice Domon et Leonie Duquet, 75205 Paris Cedex 
13, France 

2 Aalto University Metsahovi Radio Observatory, Metsahovintie 1 14, 
FIN-02540 Kylmalii, Finland 

3 African Institute for Mathematical Sciences, 6-8 Melrose Road, 
Muizenberg, Cape Town, South Africa 

4 Agenzia Spaziale Italiana Science Data Center, c/o ESRIN, via 
Galileo Galilei, Frascati, Italy 

5 Agenzia Spaziale Italiana, Viale Liegi 26, Roma, Italy 

6 Astrophysics Group, Cavendish Laboratory, University of 
Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, U.K. 

7 Astrophysics & Cosmology Research Unit, School of Mathematics, 
Statistics & Computer Science, University of KwaZulu-Natal, 
Westville Campus, Private Bag X54001, Durban 4000, South 
Africa 



8 CITA, University of Toronto, 60 St. George St., Toronto, ON M5S 
3H8, Canada 

9 CNRS, IRAP, 9 Av. colonel Roche, BP 44346, F-31028 Toulouse 
cedex 4, France 

10 California Institute of Technology, Pasadena, California, U.S.A. 

11 Centre for Theoretical Cosmology, DAMTP, University of 
Cambridge, Wilberforce Road, Cambridge CB3 0WA U.K. 

12 Centra de Estudios de Ffsica del Cosmos de Aragon (CEFCA), 
Plaza San Juan, 1, planta 2, E-44001, Teruel, Spain 

13 Computational Cosmology Center, Lawrence Berkeley National 
Laboratory, Berkeley, California, U.S.A. 

14 Consejo Superior de Investigaciones Cientfficas (CSIC), Madrid, 
Spain 

15 DSM/Irfu/SPP, CEA-Saclay, F-91 191 Gif-sur-Yvette Cedex, 
France 

16 DTU Space, National Space Institute, Technical University of 
Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby, Denmark 

17 Departement de Physique Theorique, Universite de Geneve, 24, 
Quai E. Ansermet,121 1 Geneve 4, Switzerland 

18 Departamento de Ffsica Fundamental, Facultad de Ciencias, 
Universidad de Salamanca, 37008 Salamanca, Spain 

19 Departamento de Ffsica, Universidad de Oviedo, Avda. Calvo 
Sotelo s/n, Oviedo, Spain 

20 Department of Astronomy and Astrophysics, University of 
Toronto, 50 Saint George Street, Toronto, Ontario, Canada 

21 Department of Astrophysics/IMAPP, Radboud University 
Nijmegen, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands 

22 Department of Electrical Engineering and Computer Sciences, 
University of California, Berkeley, California, U.S.A. 

23 Department of Physics & Astronomy, University of British 
Columbia, 6224 Agricultural Road, Vancouver, British Columbia, 
Canada 

24 Department of Physics and Astronomy, Dana and David Dornsife 
College of Letter, Arts and Sciences, University of Southern 
California, Los Angeles, CA 90089, U.S.A. 

25 Department of Physics and Astronomy, University College 
London, London WC IE 6BT, U.K. 

26 Department of Physics and Astronomy, University of Sussex, 
Brighton BN1 9QH, U.K. 

27 Department of Physics, Gustaf Hallstromin katu 2a, University of 
Helsinki, Helsinki, Finland 

28 Department of Physics, Princeton University, Princeton, New 
Jersey, U.S.A. 



53 



Planck Collaboration: Planck 2013 Results. XXIV. Constraints on primordial NG 



Table B.3. Summary of results from the modal estimator survey of primordial models for the main non-standard bispectrum shapes. 
This is an extended version of Table 1 1, but with results from SMICA, NILC and SEVEM. 
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Table B.4. Cross-validation of best fit feature model results for the SMICA, NILC and SEVEM foreground-cleaned maps. Results are 
only presented for feature models with better than a 2.5cr result on the full domain (see Table 12). 



Component separation . NILC SEVEM SMICA 

Feature model .... / NL ± A/ NL (cr) / NL ± A/ NL (<r) / NL ± A/ NL (cr) 



k c = 0.01125; <p = 458 ± 169 ( 2.7) 409 ± 169 ( 2.4) 434 ± 170 ( 2.6) 

k c = 0.01750; </> = -337 ± 131 (-2.6) -328 ± 128 (-2.6) -335 ± 137 (-2.4) 

k c = 0.01750; (p = 3n/4 ... 368 ± 124 ( 3.0) 348 ± 121 ( 2.9) 366 ± 126 ( 2.9) 

k c = 0.01875; <f> = -359 ± 1 18 (-3.1) -366 ± 1 15 (-3.2) -348 ± 1 18 (-3.0) 

£ c = 0.01875; = n/4 -339 ± 1 17 (-2.9) -328 ± 1 15 (-2.9) -323 ± 120 (-2.7) 

k c = 0.02000; c/> = n/4 -305 ± 1 18 (-2.6) -334 ± 118 (-2.8) -298 ± 1 19 (-2.5) 
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Table B.5. Comparison of /nl results for the hybrid polynomial and Fourier modes for a variety of non-separable and feature 
models. 



Modal basis 


Polynomial 


Fourier ISW 


Model 


/nl ± A/ N l (a - ) 


/ NL ± A/ N l (fr) 



Const (see text) 14 ± 44 ( 0.3) 10 ±44 ( 0.2) 

EFTlshape(5) 8 ± 73 ( 0.1) 6 ± 73 ( 0.1) 

EFT2shape(5) 19 + 57 ( 0.3) 13 ± 57 ( 0.2) 

DBI inflation (7) 12 + 69(0.2) -0.3 ± 70 ( 0.0) 

Ghost inflation (see text) . . -24 ± 88 (-0.3) -48 ± 89 (-0.5) 

Flat model (13) 37 ± 77 ( 0.5) 38 ± 76 ( 0.5) 

NBD (see text) 155 ± 78 ( 2.0) 116 ±92 ( 1.3) 

NBD1 flattened (14) 19 ± 13 ( 1.4) 4+19(0.2) 

NBD2 squeezed (14) 0.25 ± 0.45 ( 0.5) -0.3 ± 0.5 (-0.5) 

NBD3 non-canonical (15) . 10 ± 10 ( 1.0) 4 ± 11 ( 0.3) 

Vector model L = 1 (19) .. -5 ±47 (-0.1) -24 ± 50 (-0.5) 

VectormodelL = 2(19) .. -0.4 ±2.8 (-0.1) -1.0 ± 3.2 (-0.3) 

WarmS inflation (see text) . 1 ± 33 ( 0.04) -16 ±41 (-0.4) 

Featured = 0.015 -313 ± 144 ( 2.1) -264 ± 161 ( 1.6) 

Feature k c = 0.020 -155 ± 110 (-1.4) -167 ± 122 (-1.4) 

Feature k c = 0.025 106 ± 93 ( 1.1) 110 ±98 ( 1.1) 

Feature k c = 0.030 56 ± 89 ( 0.6) 78 ± 96 ( 0.8) 

Feature k c = 0.035 22 ± 82 ( 0.3) 15 ± 84 ( 0.2) 

Feature kc = 0.040 -0.9 ± 68 (-0.01) 4 ± 69 ( 0.1) 

Feature k c = 0.050 -0.2 ± 63 ( 0.00) 15 ± 66 ( 0.2) 

Feature k c = 0.080 16 ± 64 ( 0.2) 21±64( 0.3) 
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